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ABSTRACT 

The formation of giant planets requires accumulation of ~10 Earth mass in solids; but how do protoplanets 
acquire their mass? There are many, often competing processes that regulate the accretion rate of protoplanets. 
To assess their effects we present a new, publicly-available toy model. The rationale behind the toy model is that 
it encompasses as many physically-relevant processes as possible, but at the same time does not compromise 
its simphcity, speed, and physical insight. The toy model follows a modular structure, where key features 
- e.g., planetesimal fragmentation, radial orbital decay, nebula turbulence - can be switched on or off. Our 
model assumes three discrete components (fragments, planetesimals, and embryos) and is zero dimensional in 
space. We have tested the outcomes of the toy model against literature results and generally find satisfactory 
agreement. We include, for the first time, model features that capture the three-way interactions among small 
particles, gas, and protoplanets. Collisions among planetesimals will result in fragmentation, transferring a 
substantial amount of the solid mass to small particles, which couple strongly to the gas. Our results indicate 
that the efficiency of the accretion process then becomes very sensitive to the gas properties - especially to 
the turbulent state and the magnitude of the disk headwind (the decrease of the orbital velocity of the gas with 
respect to Keplerian) - as well as to the characteristic fragment size. 

Subject headings: planets and satellites: formation — protoplanetary disks — methods: statistical 



1. INTRODUCTION 

Giant, gas-rich planets are not formed in one day, but their 
formation should be complete within the ~10^ y r over which 

■ the g as-rich phase of protoplanetary disks last (iFedele et alJ 

. I2010h . As there are about 100 e-f oldings required to grow a 
sufficiently massive core (-10 M«: lPollacket alJIig ge!) out of 

, ISM dust grains {m ^ 10"'^ g), this nevertheless represents a 
fast process. Not only is the process fast: the planet formation 
mechanism also seems to be very efficient. Over the years sur- 

I veys like Harps and Keple r have discovered a wide number 

, and v ariety of exoplanets (iMayor et al.l 1201 It iHoward et al.l 

' 1201 ll) . Planet formation may be ubiquitous. 

The view that the initial stage of giant planet forma- 

, tion proceeds similarly to that of terrestrial planet forma- 
tion is commonly he l d. This is the core accretion paradigm 
(iMizuno et al.l Il978t iPollack et aP Il996h : a rocky core is 
formed first, whereafter it binds the nebular gas, when it be- 
comes sufficiently massive. Another paradigm for the forma- 
tion of massive planets is the disk-instability model, where 
the gaseous disk is su fficiently massive t o become gravita- 
tionally unstable (e.g., iMayeret al.ll2007t lBoleyll2009t iBossI 
I20T1I) . Here, we will consider the core accretion paradigm 
and focus exclusively on the problem of forming a ~10 M® 
core quickly enough. 

An important intermediate step is the formation of planetes- 
imals - gravitationally interacting bodi es usually thought to 
be of km-size or larger (ISafronovlll969h . These are the build- 
ing blocks for planets. Indeed, a typical assumption in pro- 
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toplanet growth studies is that the majority of solids (i.e., the 
rock and ice) of the protoplanetary disk resides in planetesi- 
mals. But how planetesimals form in the first place remains 
ill-un derstood. Models favor ing direct sticking by dust par- 
ticles (lWeidenschillingiri997l) . have recently been confronted 
with a var iety of obst acles: the fragmentation barrier at m- 
sizes (e.g..lBirnstiel et al. 201 (1); the |70uncing ba rrier at mm- 
sizes dGiittler et al EoiOtlZsoni et alj|20ial20rih ; or a charg- 
ing b arrier for even smaller (and fluffier) particles (lOkuzumil 
120091) . Perhaps these results indicate that planetesimal for- 
mation requires special conditions : e.g., sticky material prop- 
erties (like ic e: IWada et al.ll2009h; par t icle pileups in drift- 
free regions jYoudin & Chiand l2004t iKretke & LinI 120071: 
iBrauer et al.ll2008h ; concentration of chon drule-size particles 
and their subsequent gravitational collapse dCuzzi et al.ll2010l: 
Chamb ers 2010) ; or the shearing instab ility that operates on 
a den se layer of dm-m size bo ulders (lYoudin & Goodman! 
i2005tl.Tohansen et al.ll2007[l2009h . 

The transition to planetesimal size bod ies, instigates 
the runaway and ohg a rchic g rowth phases (Ida & Makinoj 
[1991 iKokubo & Idd [19961 [1991 [Ormel et al. 2010a). 
Runaway growth is fast and a protoplanet seed forms 
(IWetherill & Stewarill989h . However, at some point the run- 
away body faces a backlash from its own viscous stirring, 
thereby allowing embryo^ in neighboring zones to catch up 
in terms of mass. What follows is a two component system 
of embryos and planetesimals, with the latter containing most 

' We will employ the phi'asing 'oligarchs', 'core', 'embryo', and 'proto- 
planet' mostly as synonyms. More correctly, an embryo or protoplanet con- 
sist of a rocky (or icy) core with a gaseous atmosphere. 
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of the solid mass. The core formation challenge then amounts 
to transferring the mass reservoir of planetesimals onto these 
embryos. 

There have bee n recent reviews of the runaway/oUgarchic 
growth stag e (e.g.. G oldreich et al.ir2004l:lLevison et al.ll2010l: 
lOrmel etal.. 2010b) and we will not repeat these here, but the 
underlying problem is that during oligarchic growth embryos 
dynamically stir the planetesimals faster than they can accrete 
them. For this reason, the next e-folding in embryo mass al- 
ways takes longer than the previous. Allowing for embryo 
atmospheres (Inaba & Ikoma 2003), which enhances the ra- 
dius at which planetesimals will be captured, alleviates the 
problem somewhat, but generally modelers require massive 
disks - more massive than the minimum mass solar nebula 
(MMSN: l^denschilling|1977bl:lHavashi et al.|[T985l) bench- 
mark - to form massive protoplanets within the time that the 
gas disk disperses. 

However, planetesimal-planetesimal collisions cause a co- 
pious amount of debris, especially for massive disks. The 
question is whether these fragments - or small particles in 
general - are more conducive to growth than their km-size 
progenitors. This is not an easy question to answer since 
smaller particles interact more strongly, and in more di- 
verse ways, with the gas. On the positive side, gas damp- 
ing is much stronger, reducing their eccentricities. Simi- 
larly, the dissipative nature of gas drag enhances capture rates. 
However, orbital decay timescales can be extremely short 
jWeidenschillind ri977al) and diffusion by turbulence could 
severely dilute their midplane densities. One of the foci of this 
work is to capture these key physical processes into a common 
framework. 

The two component approximation with which oligarchic 
growth is characterized renders it ideal to be studied by 
simple numerical or even analytical models (Goldreich et alj 
2004; Chambers 2 006, 2008; Brunini & Benvenuto 200i 
Guilera et all 1201 11) . With a 'toy model' (TM) we mean 
anything where the mass distribution is approximated by 
two or three (in case fragments are included) components. 
This in contrast to models that include the full mass spec- 
trum (e.g.. ,Weiden schilling 1997; Kobavashi et al. 120101 ; 
iBromlev & Kenvon|[201 li) , which may be referred to as 'so- 
phisticated'. The goal of a TM may be to provide physi- 
cal insight or to allow for a quick (but still sufficiently ac- 
cur ate) computation. P opulation synthe s is mod els, for exam- 
ple, llda & LinI (120041) . iMordasini et~aLl (120091) . and their se- 
quels employ a TM for the embryo growth stage, if not for 
all stages. But do TMs capture all relevant physics? There is 
often a tendency for the (above) TM to include more sophisti- 
cated features, e.g., a radial dimension to follow the migration 
of protoplanets or multiple mass bins to capture the collisional 
cascade. 

In this work we will attempt to adhere more strictly to what 
we consider characterizes a TM. Specifically this demands the 
model to be; 

1. Quick. CPU speeds should be seconds, rather than 
hours, let alone days. 

2. General. That is, it should include all relevant physical 
process that could potentially affect protoplanet growth; 

3. Transparent. It must be clear which effect causes what. 

Combining the second and third properties is challenging, 
since a more complete model, characterized by more features. 



is often also a more complex one. For that reason, we have 
opted for a modular approach, where features can be switched 
on or off. As a result of the first and third property, our TM 
contains many simplification, perhaps oi'ersimplifications. 
However, in our view, a TM is complementary, not in com- 
petition, to more sophisticated models of planet formation. 
That is, the TM provides a first glance of what result may be 
expected. It can explore the parameter space most efficiently 
and assesses whether or not novel features do matter It may 
not get the details correct; its results may even be off by sev- 
eral factors of unities. But here it is where it can work in 
conjunction with sophisticated models. 

In this paper, the emphasis lies on the presentation and val- 
idation of the toy model to follow protoplanet growth. In §|2] 
we present the toy model. In §[3] we validate it against several 
literature studies. In §|4]the novel features of the toy model 
are explored. We summarize and conclude in § |5] 

2. THE TOY MODEL 

This section presents the toy model in detailQ In § |2.1| we 
give an overview of the features. In § I2.2l we describe how the 
accretion rate, or , equivalently, the embryo growth timescale, 
is computed. Sections [2. 3l and [Z4l address how the dynami- 
cal state (eccentricities of planetesimal and fragments) is ob- 
tained. Section |23] discusses collisional fragmentation, i.e., 
the way collisions affect the surface density in planetesimal 
and fragments and their characteristic sizes. Section|22] dis- 
cusses how drift motions affect the surface density of the com- 
ponents. We summarize the algorithm in § 12.71 

2.1. Overview, model switches, and caveats 

A sketch of the toy model is given in Fig.[T] The mass dis- 
tribution is approximated by three component; embryos (E), 
planetesimals (P) and fragments (F), each characterized by a 
size (sit or Re for the embryos), surface density and ec- 
centricity et where k is one of E, P, or F. These quantities are 
referred to as the state variables of the system. In our TM 
we quantify how the components influence each other. Em- 
bryos will dynamically excite the planetesimal and fragment 
population and their mass increases due to accretion. Plan- 
etesimal fragmentation transfers mass from the planetesimal 
to the fragment population. Collisions among these compo- 
nents can alter their characteristic size (sp or Sf). Finally, 
the gas from the disk plays a crucial role; it damps the ran- 
dom motions (eccentricities) of planetesimals and fragments, 
but may also excite eccentricities due to turbulence-induced 
density inhomogeneities in the gas. Furthermore, the gas de- 
pletes the surface density due to radial drift. To complement 
the model we allow (optionally) for an external inflow of mat- 
ter at a rate Mext. 

An overview of the binary switches the toy model contains 
is provided in Table[T] These can be turned on or off depen- 
dent on the problem under consideration. Due to this mod- 
ular nature, they are referred to as switches. When a switch 
equals it means the feature is turned off; when switch=l, 
it is turned on. The atmosphere enhancement (AE) and nebula 
drag (ND) switches affect the efficiency of accretion. Atmo- 
sphere enhancement results in a larger accretion cross section 
due to energy dissipation effects. Likewise, the nebular gas 
causes the orbit to be modified from the usual (dissipationless) 
2-body approximation, thereby also modifying the collision 

- The Python source code of the toy model is pubHcly available at 
http: //astro.berkeley . edu/~ormel/so£tware.html. 



Understanding how planets become massive 



3 



Turbulent 
Stirring (TS) 



Turbulent 
Stirring (TS) 



Gas drag 



^ Gas drag 




Radial drift (RD) 



External inflow 
(Mext) 

Collisions: planetesimal growth (PG) and' 
fragmentation (changes size) 

Figure 1. Sketch of processes affecting the dynamical state (red arrows) and the surface density (black arrows) of the components. Each of the thi'ee components 
is characterized by a surface density St, size, Sk (or Re for the embryo's), and excitation state (the embryos' eccentricity is by definition 0). Dashed arrows 
denote external processes like gas drag and drift motions. The random motion (eccentricity) of the planetesimal and fragment population is increased by viscous 
stining and turbulent stining, and damped by gas drag. The surface density in embryos increases due to accretion. The surface density in fragments and 
planetesimals changes to accretion, mutual collisions (erosion and fragmentation), radial drift, and inflow of solids from external regions. Collisions within the 
planetesimal and fragment population changes the characteristic size (sp and Sf) of these bodies (blue arrows). 



Table 1 

The list of binary switches (in alphabetic order) that can be turned on or off. 



Switch 


Abbr 


Atmosphere enhancement 


AE 


Erosive collisions 


ER 


Fragmentation 


FR 


Oligarchy 


OL 


Planetesimal growth 


PG 


Nebular drag 


ND 


Radial drift 


RD 


Turbulent diffusion 


TD 


Turbulent stirring 


TS 



Abbr Accounts for the . . 



Discussed in 



Increase in the capture radius due to trapping of particles in protoplanets atmospheres. 
(Fragmenting) collisions between fragments and planetesimals 
Fragmentation of planetesimals by planetesimal-planetesimal collisions 
Presence of neighboring protoplanets 

Growth of characteristic planetesimal size, due to planetesimal-planetesimal collisions 
The nebular flow in calculating the accretion rates of particles coupled to the gas. 
Depletion of the surface density due to orbital decay of planetesimals and fragments 
Vertical diffusion of particles due to turbulence, reducing their midplane density 
Additional excitation of planetesimals due to gravitational scattering by turbulence- 
induced density fluctuations 



2.2.4 



2.5.3 

IT 



2.2.1 



2.5.4 



2.2.3 



andlTel 



2.2.3 



2.3.2 



rate. Fragmentation among planetesimals is included when 
FR=1, and planetesimal-fragment erosive collisions are addi- 
tionally included when ER=1. The oligarchic switch (OL) indi- 
cates the presence of neighboring protoplanets, which affects 
the accretion rate and the drift timescales of solids. Growth 
of planetesimals is taken into account when PG=1; otherwise 
the planetesimals characteristic size sp does not evolve. Ra- 
dial drift (by fragments and planetesimals) is included when 
RD=1. Stirring due to turbulence-induced density fluctuations 
is included when TS=1; otherwise, viscous stirring by em- 
bryos is the only mechanism that excites the eccentricities of 
the planetesimal bodies. Finally, we consider the possibil- 
ity of vertical diffusion of particles due to turbulent diffusion. 
The TD switch is then turned on. 

Other (continuous) parameters of the toy model are listed 
in Table|2] The toy model is further characterized by the fol- 
lowing properties and assumptions: 

1 . We prefer physically-motivated order of magnitude ex- 
pressions rather than detailed black-box formulas. This 
implies that we favor physical insights at the expense of 
precision. 

2. There is no explicit radial dimension. Type-I migration 
is neglected. Removal of (small) particles by aerody- 
namic drift is included (if RD=1), and an (ad-hoc) in- 
jection rate of mass may be pr ovided (Mext 0). A 
global interpretation (see § 12.61 ) can be made when oli- 
garchy is assumed (0L=1), but fundamentally the model 



is strictly local. 

The model assumes a bimodal (or trimodal) distribution 
for the mass. This implies an instantaneous jump from 
the planetesimal size sp to the fragment size Sf. Fur- 
thermore, the three-component assumption is not appli- 
cable to follow the runaway growth phase (the phase 
preceding oligarchic growth). 

Regarding the dynamical state (computation of eccen- 
tricities), the model assumes a balance between the 
(viscous) stirring timescale and a limiting timescale 
(§ 12.41 ). We do not integrate the stirring rates (deceit or 
dif/dt). This setup implies memoryless behavior; e.g., 
the eccentricities follow from the current values of 
and Me- 

We implicitly assume that <sc -H I,p, such that the 
random motions of protoplanets are negligible due to 
dynamical friction (e^ - 0). The toy model is therefore 
only valid as long as oUgarchy pertains. 



2.2. The accretion rate 
2.2.1. Hill units and key definitions 
We define the total accretion rate of the embryo as 

dME 
dt 



(1) 



/t=P,F 
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Table 2 

List of pai'ameters and frequently used symbols. 



Symbol" Default value'' Description 



Reference 







Relative velocity between particles of component k and / 






Gas surface density CLg); solid surface density of component k (Sj) 


^ini 


(mmsn) 


Initial total suiface density in solids 


n 


(mmsn) 


Local orbital frequency corresponding to disk radius ciq 


as 




Embryo radius in terms of Hill radius 


0!ss 


10-"* 


IShakura & SunvaevI (19731) turbulence viscositv parameter 


7 


1.4 


Adiabatic constant of the gas 




0.1 


Control parameter for timestep increment 


K 


0.01 cm^ g"' 


Atmosphere opacity 


P 




Gas density within the embryo's atmosphere (Appendix) 


Plc.s] 


1 g cm"^ 


Internal density of the embryo's core (pc); and solids (planetesimals, fragments; p, ) 


Pg 


(mmsn) 


Gas density of the nebula midplane at a disk radius ag 


Ckl 




Cross section for collisions between components k and / (main paper); dimensionless 
density (Appendix) 


Tfr 


(mmsn) 


Dimensionless friction time (= Tf,-f2) 


Cacc 


1 or 1.5 


Accretion rate correction factor due to embryo-embryo merging 


Cdrift 


0.5 


Prefactor for determining the drift timescale, applicable when 0L= 1 


Ccol 


36 


Numerical constant that enters the stiiTing rate expression 


^gg 


9.0 


Prefactor in strength expression for gravitationally-bound bodies 


^9 




Scaleheight gas disk (= c , / Q.) 




1 Lo 


Stellar luminosity 




1 Mo 


Stellar mass 






[core. Embryo] mass"^ 




10-* Me 


Initial core mass 






Isolation mass 







Mass accretion rate of solids from exterior regions drifting past the embryo 






Specific collision rate in terms of Hill units 






Specific collision rate in Hill units in the 2-dimensional limit 


^'vs 




Viscous stirring rate in Hill units 


^vs,sd 


73 


Max stirring rate, obtained at Hill eccentricities e;, <s: 1 


R[a,b.c,h.El 




[Atmosphere-enhanced, Bondi, core. Hill, embryo] radius"^ 


Tcxil.kk 




Collision timescale for planetesimal-planetesimal or fragment-fragment collisions 


^[E,dr,tr],<: 




[Depletion, drift, friction] timescale of planetesimals or fragments 






Growth timescale of embryos due to accretion of component k 




(mmsn) 


Temperature nebula at a = ciq 


T'stir 




Net stirring timescale (maximum of Tys and Tis) 


7'[ts,vs] 

Gog 




Timescale for [turbulent, viscous] stiiTing 


10' cm^ s-^ 


Prefactor for the gravity regime in the Q'^ law 


Gos 

Grf 


2.1 cm' s^^g-' 


Prefactor for the strength regime in the law 
Material strength 


w. 




Radial width of the simulated region 






Dimensionless parameter for the calculation of atmosphere structure 




5.0 AU 


Disk radius (semi-major axis) 


b 


12.5 


Separation of oligarchs in units of the Hill radius Ri, 




[1.19, -0.45] 


Exponent of [gravity, strength] term in Q'^ law 


c., 


(mmsn) 


Isothermal sound speed 






[Orbital, Hill eccentricity] of component k 


eiii 


2.34 


Transition Hill eccentricity for Pys fit 


hk 




Scaleheight particle layer accounting for Keplerian inclination and turbulent dilfusion 


k 




Orbital inclination 


Ik 




Ratio collision energy over collision strength 


Sk 




Radius planetesimals or fragments 


^F.ini 




Initial radius fragments 


*P,ini 


10 km 


Initial radius planetesimals 


VK 




Keplerian (orbital) velocity corresponding to oq 


Vesc 




Escape velocity of the embryo 


Vh 




Hill velocity 


"hw 


54 m s 


Headwind velocity experienced a by heavy body due to sub-Keplerian moving gas 


t>td,k 




Radial drift velocity for particles of component k 



2.2.3 



Section l2.2.2 
Equation )49 
Equatio n (49 
Section l2.Zl 
Equation (3) 
Section 
Section 
Section 
Section l233l 
Appendix IB I 



Section [3] 

§ [2T2l AppendixH] 



2.2.4 



Section sl2.2.3l and l2.3.3l 
Section [2.2. II 
Section |2. 6 1 equation (44) 
Equation (5) 
Equation j29) 

Section[3] 

Equation j50t 



Section[3] 

Equation j5H 



2.2.1 



2.2.3 



Section 
Section 
Section 
Equatio n fl5 

Section l2.2.3 

Sections l2Tu and l2.2.4l 

Equation i26\ 

Equations'^, (19), and g2) 
Equation 1 4} 
Equation 1 50) 
Equation fl8 
Equations Tl7t 
Equation i29\ 
Equation (29) 
Equation (29) 
Section l2.6l 
Equation (13t 



i and fie) 



2.2.1 



Section jlXT 
Equation (29 
Section 3[ 
Section ~ 
Section 
Equatio n j9 
Section l2.2.2l 
Equation (28| 



2!T7l Figure[3] 



Figure |5] 



2.2.2 



2.2.1 



Section 

Section 

Section |2.3.3. 
Equations (2 It 



Equation (lit 
■ and<22) 



" The dummy index k (one of E, P, or F) indicates one of the three components of the toy model (embryos, planetesimals, or fragments). Occasionally, 
similar symbols are combined in one fine using squai'e brackets. 

A number indicates the default value of the parameter; '(mmsn)' implies that (by default) we follow the minimum-mass solar nebula model (§[3). 

^ The atmosphere is assumed not to contribute significantly to the embryo's mass, i.e.. Me = M^. The embryo radius Re, as it appears throughout S I2.2l is 
either the atmosphere capture radius So (when AE=1) or the core radius R^ (when AE=S). 
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where Me is the embryo mass, Cacc a correction factor for 
oligarchic growth (discussed below), the surface density 
and the sum is over the planetesimals (k - P) and fragments 
(k - F) components. The Hill radius Ri, is defined as 



Rh 



ao 



Me 
3M, 



1/3 



(2) 



with flo the semi-major axis of the protoplanet, M^, the stellar 
mass, ui, = Ri,Q. the Hill velocity, and Q. = Q(fl()) the orbital 
(Kepler) frequency at oq. Finally, Pcoi is the specific accretion 
rate, which depends on the eccentricities (e) of the particles, 
their scaleheight (or inclination, /), and the dimensionless ra- 
dius of the embryo: 



Ue 



Re 



1 /9M, 



flo \4-7Tpc 



1/3 



10" 



5 AU 



3 g cm 



-1/3 



(3) 

where Re is the radius of the embryo (the core radius) and pc 
its density. Note that oe is primarily a function of disk radius 
flQ. We further introduce the Hill eccentricity e/, = ey^/y;,. 
The use of the Hill eccentricity e/, is convenient since it dis- 
tinguishes the dispersion-dominated (d.d.) regime (e;, > 1), 
where approach velocities are determined by the planetesi- 
mal eccentricities, from the shear-dominated (s.d.) regime 
(e;, < 1), where approach velocities are given by the shear 
of the Keplerian-rotating disk. 

In equation ([T]) Cacc represents a correction factor to dM/dt, 
relevant for oligarchic growth, due to the self-accretion of em- 
bryos. Its value can be obtained as follows. During their 
rowth embryos keep a distance of bRi, with h ~ 12.5 constant 
Kokubo & Ida 1998 ). Since Ri, cc M^'^ the number of em- 
bryos has halved when their mass has increased by a factor 8. 
Thus, the growth by a factor of 8 is split with embryo-embryo 
growth accounting for a factor of 2 and embryo-planetesimal 
growth accounting for a factor of 4. Therefore, Cacc = 1-5, 
when we consider oligarchic growth (0L=1). Otherwise, if we 
consider a single protoplanet (OL=0) the correction does not 
apply and Cacc = 1- 

From equation ([T]i, the accretion (or growth) timescale for 
each component is defined as 



Ter.k = 



Me 



dME/dt C^ccPco\,ARhVh 



(4) 



Next, we derive expressions for f coi for the d.d.- (e;, > 1) and 
s.d. -regimes (e/, < 1) . 

2.2.2. The dispersion-dominated regime: e;, > 1 

In the dispersion-dominated (d.d.) regime the relative ve- 
locity between particle and protoplanet is given by the eccen- 
tricity of the former, A^ep = epVK where vk - aoQ. is the 
Kepler velocity. The d.d. -regime is generally valid for big 
bodies, orbiting on a well defined Kepler orbit of inclination 
that is half of the eccentricity, / - e/2- the so called equilib- 
rium value. For Pcoi we adopt 



~ ^col -» 



(5) 



with Ccoi = 36 (iGreenzweig & Lissaueii[T990l[T992h . 

As an illustration, let us verify the validity of equation ^ 
by an order-of-magnitude argument. The rate at which a pro- 
toplanet sweeps up particles can be written as 

dME 

—— a npo-EpiAvEp)mp, (6) 
dt 



where («/)crEpAL>Ep)"' is the collision timescale between plan- 
etesimals and embryos, nip the mass of the planetesimal, np 
the number density of planetesimals, ctep the gravitationally- 
enlarged collision cross section. Gravitational focusing en- 
larges the collisional cross section to ctep - irR^iugsc/ ^^ep)' 
where Hesc - ^I2GMe I Re is the escape velocity of the em- 
bryo with G Newton's gravitational constantQ Furthermore, 
we can write np - I,p/2hpmp, with hp x iciQ the scaleheight 
of the planetesimal layer and / - e/2 the inclination. Then, 
equation ^ becomes 



dME 
"dT 



r>2 \ ^\ I '-^esc 



- 

evK 



(7) 



It can be shown that RE^i^c - (^Rhv\\ and since e;,y/, = cvk 
the above equation, up to a constant, equals equation ([T) with 
equation (|5]l for Pcoi- The larger value of Ccoi above is due to 
the averaging over a distribution of eccentricities. 

2.2.3. The shear-dominated regime: ei, < 1 

In the shear-dominated (s.d.) regime - valid mostly for 
smaller particles - the expression for Pcoi is not so straight- 
forward as equation First, inclina tions start to decouple 
from the eccentricities jlda et al.lll993h : the ratio //e = 1/2 is 
no longer maintained. Nevertheless, we will assume this rela- 
tion to hold as we do not separately compute the inclinations. 
More important is the way the gas interacts with the (small) 
fragments - a subtlety that is sometimes neglected in studies 
of oligarchic growth. We consider two effects: 

1 . When disks are turbulent, diffusion of particles causes 
them to spread over a vertical height that may be 
much larger than the corresponding value due to their 
inclination (/flo). Here, we take the ensuing scaleheight 
for the dust as 



Ha 



O's: 



(8) 



(lYoudin & Lithwickll2007h where Tf, = TfrQ is the di - 
mensionless stopping time (see definition in § 12.3.3b . 
ass the fa miliar turbulent diffusion pa rameter first intro- 
duced bv lShakura & SunvaevI ( 119731) . and Hg the scale- 
height of the gas. The net scaleheight of the particle 
layer - the quantity used in Pcoi (see below) - is ap- 
proximated as: 



hp - max /flo 



[''^9 \l ~ — 1 • 

■' V ttss +Tfr/ 



(TD=1) 



(9) 



Thus, when Tfr <^ o-ss the particles' scaleheight be- 
comes equal to that of the gas. But even when Tf,- > a^s 
turbulent diffusion can result in a scaleheight much 
larger than that from inclination alone {hp - iao when 
Q^ss = or TD=0). When gas-drag effects do not affect 
Pco\ (i.e., ND=0), the maximum collision rate is attained 
when the particles reside in a thin layer {hp 0). For 

1 /'^ 

this g as free, 2D configuratio n Pcoi — ^'coi,2d ~ Ho'g 
(e.g., Ilda & Nakazawalll989l) . This 2D limit is vaHd 
when the particle scaleheight is less than bco\, the im- 
pact radius for accretion. Otherwise, the local density 

' More coiTectly, the focusing factor is 1 + vI^^/{Avep)^ but in our case 
we ignore the unity term since gravitational focusing factors for ohgarchic 
growth are always 2> 1 . 
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Figure 2. Specific collision rate Pcol for a 0. 1 protoplanet situated at 5 AU. (left) P^a] as function of Hill eccentricity, e;, . For a gas-free system (ND=TD=S) 
Pcol is a function of e;, only (thick solid line). Adding turbulent diffusion (TD=1 with cj-js = IQ-^; dashed lines) lowers the collision rate in the shear-dominated 
regime and renders it a function of the particle size or dimensionless stopping time Tfp = Tf^Cl. rather than of e;, . ( right) Pcol as function of Tfj in the s.d. -regime, 
valid for <*; 1 and a headwind of «hw = 54ms-'. The curves give P^al for the case including turbulent diffusion (dashed curves) and for the modification of 
the collisional cross section by nebular drag (gray curves; ND=1). 



of particles is diluted by a factor bco\/hp and the colli- 
sion rate reduces accordingly: 



Pco\ = ^"€01,2(1 min 1 1 , ^ j , 



(10) 



1/2 

where in the s.d . -regime, br.„] - 1-7q '£- Rh (see e.g., 
Ilnaba et al1l200Tt loimel & Klahill20Tol) . 

2. Gas drag forces cause an encounter to deviate from 
its 2-body, energy-conserving, trajectory. Consider- 
ing this effect (ND=1) gives rise to a Pcoi and ^coi that 
are different from the gas-free case, discussed above. 
Ormel & KlahJ (120101) (see also lPerets & Murrav-Clavl 



201 Ih have accounted for the behavior of the gas drag 



during the encounter and have derived analytic expres- 
sions for bco\ and P^oi in the presence of gas drag. 
Therefore, when the nebular drag swit ch is turned on 
ND=1), ^coi and Pcoi,2d follow from the lOrmel & Klahil 
1201 Oh study. These expressions are summarized in 
Appendix IaI If, in addition, vertical diffusion is con- 
sidered (TD=1), the correction factor to Pcoi,2d (eq- ifTOl ) 
is also applied. 



In Fig. |2] the effects of these features on Pcoi is illustrated 
for a 0.1 Me protoplanet at 5 AU of internal density pc - 
3 g cm-"*. The thick solid line in Fig. [2^ shows Pcoi for the 
gas-free case, which is equivalent to setting the switches ND 
and TD to zero. Pcoi increases with decreasing eccentricity. In 
the d.d. -regime (e/, > 1) it is given by equation (|5]l. The s.d.- 
regime splits into two regimes: (i) the 2D limit corresponding 
to a thin particle layer (e/, «: 1; dotted horizontal line); and 
a transition regime, where the thickness of the particle layer 
is larger than the accretion radius, iao > bco\- (Note that we 
assume i = e/2 throughout). Although our formulation is 
slightly different, the Pcoi we obtain for the gas-free limit is 
consistent with the expressions given by Inaba et al.. (,2001.). 

Allowing (small) particles to diffuse by turbulence (TD=1) 
decreases their density at the midplane and decreases Pcoi- 
The dashed lines in Fig.|2^ show that the corresponding re- 
duction, valid for a^^ - 10 -*, can become quite significant for 
small particles. The correction factor is only applied in the 



s.d. -regime, i.e., for e/, < 1, since it only matters for small par- 
ticles that anyway satisfy e/, <*: 1 . For this reason, it is more 
useful to plot Pcoi as function of size, or, rather, dimensionless 
friction time, see Fig.|2j5, black dashed line. With increasing 
Tfr, Pcoi increases until the 2D-limit is reached (dotted hori- 
zontal line). 

Accounting for the nebular flow around the protoplanet 
(ND=1) results in accretion rates that are quit e different from 
Pcni fo r the gas free case. As outlined by lOrmel & Klahrl 
( 120101) an important parameter in determining the accretion 
rates is the magnitude of the disk headwind Uhw The gas in the 
protoplanetary disk is partly pressure supported, resulting in a 
angular velocity slow er than the gas by a quantity = J]Vk, 
where tj is defined as (lAdachi et al.lll976h : 



Tjiao) 



1 



dP„ 



IpgOoD? da 



I El 

'2[vl 



dlnPg 
51nfl 



(11) 



where pg is the gas density of the nebula at disk radius ao, 
Pg — Pgc] the gas pressure, the isothermal sound speed, and 
the gradient is evaluated at ao. A heavy body (planet or plan- 
etesimal) is not influenced by the pressure gradient and moves 
at the Keplerian velocity; it therefore experiences a headwind 
of magnitude - t]Vk- The solid gray line in Fig.[2}5 gives 
Pco\ for a headwind of Uhw - 54 m s-' (see AppendixlAlfor a 
summary how Pcoi is derived when ND=1). For Tf, » 1 one re- 
covers the 2D, gas-free limit (large bodies are not influenced 
by gas drag). Somewhat smaller particles (xfr > 1), however, 
have a finite probability to be captured during their passage 
through the Hill sphere. This capture mechanism is indepen- 
dent of the physical size of the protoplanet. If the protoplanet 
is sufficiently massive, it also operates for particles Tfr < 1 . 
However, for very small particles (rfr <K 1) the gas drag force 
is so strong that accretion is suppressed, because particles are 
virtually glued to the gas. In that case, the accretion rate falls 
below that of the gas free-limit. Including turbulent diffusion 
(dashed gray line) further decreases Pcoi- 

2.2.4. Atmosphere enhancement (AE) of the collision radii 

The collision radii further increase when protoplanets ac- 
quire an atmosphere. Since this is generally true, i.e., for 
fragments as well as planetesimals, we discuss it separately. 
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However, its effects are most pronounced for small particles. 

A large protoplanet will acquire a dense atmosphere, 
which enhance s the e ffect ive capture radius for colli sions. 
Ilnaba & Ikomal (120031) and iTanigawa & Ohtsukil (120101) have 
studied this effect by calculating the energy loss of the plan- 
etesimals near the point of their closest approach to the planet 
(where the atmosphere is densest). In this way, they derive an 
enhanced radius for accretion, Rg. To obtain the correct ac- 
cretion rate, this 'atmosphere radius' must be used instead of 
the core radius {Rc or the dimensionless ac) in all of the above 
expressions. This is the procedure we adopt here. In any case, 
we cap Ra when it reaches either the Hill radius or the Bondi 
ra dius, rI - GM/ycl wi th y = 1.4. 

Ilnaba & Ikomal (l2003h provide a solution for the density 
structure of the atmosphere. Their solution is quite general; 
it allows for multiple shells, where the conditions regard- 
ing, e.g., the opacity can be different. Here, we will make 
the additional simplification that the opacity is low and con- 
stant, K - 0.01 cm g"', reflecting essentially a grain-free at- 
mosphere. A low opacity may find some justification from 
the fact that (i) most of the mass is in macroscopic bodies 
(planetesimals); and (ii) grains may otherwise quickly coag- 
ulate (Movshovitz et al. 2010; Helled & Bodenheimer 201 1|. 
Although, a constant opacity throughout the atmosphere is 
arguably a somewhat crude assumption, our simplified ap- 
proach man ages to obtain good co rrespondence with the full 
solutions of Ilnaba & Ikomal (l2003h . see Appendix IbI 

By virtue of the low opacity, we assume that the atmosphere 
is always radiative. In Appendix|B]we approximate the equa- 
tions of stellar structure and obtain an explicit, closed-form 
solution for the density as function of depth, p = p(r). Us- 
ing this solution we straightforwardly solve for the enhanced 
radius, 7?„: 



R„ a! Rh 



1 



2Wneh(o-a - 1) + log Cr„ 



1 

+ 

Xl 



7 

4(4Wneb)'^' 

r 



; (1 < cr„ < o-i) 
(12) 

where cr,, = Pa/Pg is the ratio of the density required to cap- 
ture the particle and the nebula gas density (see AppendixlBli. 
o"! = 1 /5Wneb a dimensionless density threshold that specifies 
the transition from the nearly isothermal regime to the pres- 
sure supported regime, and xi = R\/Rb = 1 + 2Wneb(c'"i - 
1 ) + log (Ti the dimensionless radius corresponding to Ri. The 
parameter that sets the structure of the atmosphere is Wneb- 



3KLr 



ksKpgPcRl 
l6o-^bl^T^jTgi- 



UO^km/ 
Pc Pg 



10 ^ cm^ g ' 



(13) 



(14) 



gem 



10-10 g cm-3 



100 K 



-3 



\Myr 



where kg is Boltzmann's constant, p. - 2.34 the mean molec- 
ular mass, cTsi, Stefan-Boltzmann constant, Pg - kTgp/p the 
pressure of the nebula, and Tg the nebular temperature. For 
the luminosity due to planetesimal accretion L^. we have sub- 
stituted Lc = {GMJRc)dMJdt = GMIIRJ^,, with T^, the 
growth timescale and Rc the core radius. Our solution for Ra 
(eq. |[T2l ) requires that Wneb ^ 1- 
Admittedly, a certain level of ambiguity is present between 
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Figure 3. Comparison between the adopted stilting curve of our toy model 
(dashed line) to the fit of Ohtsuki et al. (2002) {solid curve). A 1:2 ratio be- 
tween inclinations and eccentricities is assumed. Our two-segment approxi- 
mation breaks at an eccentricity e;, = e;,i = 2.34. 

the atmosphere-enhanced radius Rg and the enhanced colli- 
sion radius, /^coIj 

obtained from lOrmel & Klahii ( 120101) , due 
to nebular gas-drag effects. Both s tudies account for gas drag; 
but where lOrmel & Klahri ( 120101) consider the change (en- 
hancement or decrease) of f coi to originate from the inter- 
action w ith the nebular ^ as flow (assumed to have a constant 
densitv). Ilnaba & Ikomal ( l2003h calculate the enhancement of 
Pcoi due to the increase in density in the direct vicinity of 
the planet, i.e., its atmosphere, while ignoring the flow pat- 
tern. Indeed, particles strongly coupled to the gas flow may 
be inhibited to accrete (Se kiva & Takedall2003l; IPaardekoopeii 
,2007; Ormel & Klahr 201^ 

Further work must account si- 
multaneously for both the gas flow and the density structure 
of the planet - like in iLissauer et alJ (l2009h . but for a wide 
parameter range. 

In this study, we will merge the Ilnaba & Ikomal (120031) and 
lOrmel & Klahr (2010) approaches. Specifically, we first cal- 
culate the atmosphere radius Rg (eq. lfT2l ). which then re- 
places the embryo radius Re in all the above expressions in- 
volving f col- That means that equation (O is superseded by 
qe = Ra/Rh- For the protoplanetary growth stage, the atmo- 
sphere does not contribute much to the total mass; i.e., the 
embryo mass Me remains dominated by the mass of the solid 
core Mc- 

2.3. Stirring, friction, and drift timescales 
2.3.1. Viscous stirring 

Similar to the accretion rate, equation ([TJ, we define the stir- 
ring rate for (planetesimal) bodies also in terms of Hill units: 



de^ 

-A = P,,NERhUh 
at 



(15) 



where Ne is the column density of embryos and Pvs the spe- 
cific viscous stirring rate, which is function of eccentricity 
(ei,). Figurell] shows th e dimensionless viscous stirring rate 
Pvs dOhtsuki et al.ll2002l) for the case that inclinations and ec- 
centricities hav e reached a so called equilibrium [i/e x 0.5; 
llda et al.lll993l) . The basic features - a plateau for e;, « 1 and 
a e^-^ decline fo r g/? » 1 - can be u nderstood from geomet- 
rica'l ar guments (IIdalll990l;[Goldreich et al. 2004; Ormel et all 
l2010bh . but the detailed functional form of Pyajeh) follows 
from calibration against A^ -body experiments (IStewart & Idal 
120001; lOhtsuki et all l2002h . Here, we simpli fy matters by 
adopting a simple 2-segment fit-by-eye to the lOhtsuki et alj 
(l2002i) curve (which is itself a fit), accurate enough for our 
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Figure 4. Viscous {black lines) and turbulent (gray lines) stirring timescales 
for protoplanets of 0. 1 and 10"' at a distance of 5AU. 

level of precision. See the dashed line in Fig. [3] 
Using equation (fTSI ), the viscous stirring timescale is 



deitldt 



dejjdt 



Anhelao 



(16) 



where we have inserted Ne - l/ilnaobRi,). That is, each 
protoplanet is assumed to stir an annulus of width b times its 
Hill radius of the protoplanet, where S is a parameter whose 
defauh value is ^ = 12.5 jKokubo & Idalll998h . The stirring 
timescale is plotted by black lines in Fig. |4]for conditions typ- 
ical of 5 AU. 

2.3.2. Turbulent stirring 

Apart from stirring by protoplanets, we optionally in- 
clude stirring of planetesimals by turbulent-induced density 
fluctuations in the gas. This 'turbulent stir ring' (TS) has 
recent l y been explored by (among othe rs) Laughl in et al. 
TdaetaLl d2008l): iBaruteau & LinI ( I201Q) : .Yang et al. 
Gressel et al.l (1201 ll) . We follow here a prescription 



(2004); 



2009h: 



bv IBaruteau & LinI (120101) . which connects the intensity of 
the gravitational torques that the fluctuations induce to the 
turbulent o-ss parameter. This results in a turbulent stirring 
timescale of (see Appendix ICli: 



5 X lO-^e^ /// floS, 



(17) 



Turbulent stirring is most important for small protoplanets and 
massive turbulent disks (see Fig.|4]i. Otherwise, viscous stir- 
ring dominates (see Appendix ICl and Fig.|4]i. 

When turbulent stirring is implemented, we compute both 
Jvs and Tts. Otherwise, we only use viscous stirring, i.e.. 



min(rv 



(TS=Q) 

:;rts) (TS=i) 



(18) 



2.3.3. Gas friction and radial drift 



The stopping or friction time - the time over which the par- 
ticle velocity is changed by friction with the gas - is defined 
as: 



rfr,A-(y!;as) 



3CdPgVgns,k 
VthPg 



(Sk < I'mfp); 



(19) 



(k - P, F) where is the radius of the particle, ps its inter- 
nal density, /mfp the mean free path of the gas, Cd the drag 



constant, which is a function of particle size, Ugas.A- the gas- 
particle relative velocit y, and Vth the thermal spe ed. For the 
drag constant we follow lWeidenschillingl ( Il977al) and use 



r - 



0.44 



if Rep < 1 

if l<Rep< 800 

if Rep > 800 



(20) 



with Rep = 



2sv„ 



gas/i'moi the particle Reynolds number with 
Vmo\ - Wpi'th/2 the molecular viscosity of the gas. 

Due to the disk headwind fhw (see eq. ITTl ). solids will lose 
angular momentum. This causes a particle to spiral in at a 
radial velocity (e.g .IWeidenschiUinglll977al : lN"akagawa et alJ 
[T986tlBraueretalJl2007h : 



2Tf, 



l+r; 



(21) 



The above expression holds for small particles: they will pile 
up in pressure maximums where Hhw = 0. Equation ( ISTT i is 
zero if the headwind vanishes - but this does not fully capture 
the behavior of heavy bodies on eccentric orbits. To first or- 
der in eccentricity (the epicyclic approximation), they feel a 
headwind 50% of the time and a tailwind for the other 50%, 
such that the net torque the gas drag exerts is zero. However, 
to second order in e there is a net effect: angular momentum 
is removed from the body, causing it to drift in. This effect 
amounts to the planetesimal experiencing a net headwind of 
~e~VK and we write, more generally than equation (l2Tl i 



(22) 



(Within order of unity, this simple formula ca ptures the de- 
tailed orbit-averaging calculations performed bv lAdachi et al.l 

CM) 

The particle-gas velocity Ugas is constructed in a similar 
vein. It consist of a part due to the inclinations and eccentric- 
ities of the bodies (^evx) and a contribution from the slower 
than Keplerian moving gas. The latter is due to the radial drift 
(eq. ETI ) when Tfi 1 and equals (the azimuthal head- 
wind) when Tfr » 1 . Therefore we have approximately 



2rfr + t\ 
' 1+t2 ■ 

rr 



(23) 



The friction time can be found using an iterative procedure 
involving equations (fT9b and ( 1231 ). 

2.4. Solving for e/, 

The adopted power-law approximations for the stirring 
timescale readily allow us to solve for the Hill eccentric- 
ity assuming that the stirring is balanced by another mech- 
anism. This determines our two growth modes: (i) in equilib- 
rium growth Tstir is balanced by the friction time r^, whereas 
(ii) in non-equilibrium growth Tstii is balanced by the growth 
timescale of the protoplanet, Tgiow 

2.4.1. Equilibrium growth: T^i„ = Tf,. 

In this regime Tstii = Tf, < Tg,. In the typical case that 
both Tsti, and Tfr depend on eccentricity, we use an iterative 
procedure to solve for e/,. For example: 

0. Start with an initial value of Ugas, e.g., Ugas = fhw • 

1. Usin g Lig as, obtain the friction time according to equa- 
tion (fTgf. 



Understanding how planets become massive 



9 



2. Equate Tsnr - Tf^ and solve for the eccentricity e/, by 
means of the two segment fit presented in Fig.|4] 

3. Determine a new Ugas using equation i23[ and the value 
of e/, derived in the previous step. 

We iterate the loop until convergence is achieved. We then 
obtain e;,,A and Tgi-j^ via equation (|4]i. Finally, we check for 
the validity of the equilibrium regime: of the three timescales 
involved, T^^k must be the longest. If this condition is not 
satisfied, we turn to the nonequilibrium mode instead. 

2.4.2. Nonequilibrium growth: T^nr = Tg,. 
In this regime Tstii - Tgi- < Tf^. We numerically solve 



(24) 



for e/,. We then obtain the (Hill) eccentricity consistent with 
the non-equilibrium regime. For simplicity, when solving 
the above equation we assume that Tgi- is independent of the 
friction time (i.e., we ignore drag or atmosphere enhance- 
ment effects). This is justified, since the validity of the non- 
equilibrium solution is limited to the initial stages (low T^) 
and big planetesimals (large Tfr). 

2.4.3. Friction regime: Tf, = 

When both the equilibrium and the nonequilibrium regime 
are inapplicable, we obtain e/, consistent with the friction 
regime: rf, = Tgi < Tstii. We found that the friction regime 
may become important after planetesimals had been excited 
to a large eccentricity by turbulent stirring. 

2.5. Collisional fragmentation 
2.5.1. Material strength and the collision timescale 

The collision timescale for particles within the same com- 
ponent is given by „ = nkCTkki^Vkk (k = P, F), where «^ is 
the number density, (Jkk the collision cross section, and Ay^-A 
the relative velocity between two ^-particles. We relate «i to 
Y.k via 



2hk 



with hk the scaleheight of the particle layer (eq. ||9l)- 
equal-size spheres, crkk - 47ri^ and 



! cohkk - 



2psSkhk 

SI-kAvkk ' 



(25) 
For 

(26) 



For 'heavy' particles like planetesimals the relative velocity 
and scaleheight are determined by the Keplerian elements ek 
and ik- Then, Avkk = SkUK ~ likaoQ = 2/1^0 and 



col.K 



(Tfr » 1). 



(27) 



In the toy model we assume that the fragmentation rate is 
determined by the ratio of the specific collision energy and the 
specific strength of the material: 



1 

qk = -mf^iAvkkY 



(mi + m2)Qj 



(28) 



d,k 



where m\,m2 are the masses of the collision partners and m^ 
the reduced mass. For equal-mass particles m\ - m2 = 2ot^ - 




10' 10"* 10' 
Particle size [cm] 

Figure 5. Specific strength of bodies as function of their size (solid curve). 
The material properties reflect those of ice. The thin dashed, dashed-dotted, 
and dotted lines denote the three contributions out of which Qj is composed. 
The arrow indicates the initial fragment size. 



mk- The material strength Q*^ is assumed to obey the following 
size dependence: 

e:(^) - Go. (j^f + Qo,P. (j^f + C,,vUs), (29) 

where the three terms on the RHS respectively denote con- 
tributions from the strength regime, the gravity regime, and 
the gravitational potential. The constant Cgg is fixed at 9 
dStewart & Leinhardtll2009l) . I n Fig.[5]we plotg as fu nction 
of size for ice-like materials jBenz & AsphaudU999l) . The 
corresponding values are also listed in Table |2l 

The fragmentation among planetesimals of size sp in- 
creases the abundance of fragments of size Sf <K sp. The 
intermediate scales - i.e., the collisional cascade - is not in- 
cluded in the toy model. Erosive collisions between fragments 
and planetesimals can be included (ER=1), which results in an 
increase of the fragmentation rate. We now provide expres- 
sions for the fragmentation rate ipf due to mutual collisions 
among planetesimals and erosive collisions with fragments, 
respectively. 

2.5.2. Planetesimal-planetesimal collisions 

We adopt the prescription of iKobayashi & Tanakal (120101) 
for the excavated mass (m^) in a fragmenting collision: 



1 



-(mi + mo). 



where A is defined as 



{mf,(Avf 
(mi + m2)Qd' 



(30) 



(31) 



For equal-mass planetesimals, cp - (pp - (Aypp)^/82d,p = qp- 
Collisions take place on a planetesimal collision timescale 
7"coi,pp (eq- 1271 ). during which each body loses a fraction 
qpl(l + qp) of its mass, resulting in a fragmentation rate of 



dY.p 



qp 



(1 + qp) Tcoi,?? 



(32) 



In the limit of qp » 1 equation (l32l i reduces io dl.pl dt - 
-Sp/Tcoi,?? and the surface density in planetesimals is re- 
duced on a collision timescale. For qp ^ 1 the depletion 
timescale increases by a factor . 
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2.5.3. Planetesimal-fragment collisions (erosion) 

For planetesimal-fragment collisions we can assume that 
the planetesimals determine the relative velocity (Aypp - 
Aiipp) and cross section, crpp = ns~p. Equation ( [3T1 ) then reads 



« 1. 



(33) 



The relevant collision timescale is now the timescale on which 
a fragment collides with a planetesimal 

■r-I 



= nefF(A(;pF)crpF. 



(34) 



The above equation must account for the fact that the two pop- 
ulations have different scaleheights: if hf » hp (resulting, 
for example, from turbulent diffusion) the density is set by the 
scaleheight of the fragments; otherwise (if hf <^ hp) by that 
of the planetesimals. Thus, «efF, the effective density of frag- 
ments for the collision, involves the largest scaleheight and 
we write yi^s - /2/!niax'«f with h^ia - ina\{hf, hp): 



- col.PF 



(35) 



where we used Awpp - IhpD.. During the collision, the plan- 
etesimal loses a fraction ^ (eq. ||33l ) of its mass. The frag- 
mentation rate then reads 



dLp 
dt 



Aqpmp'Lpns^pI.F^ I h 



2"col,PF 

qpl.f I hf 



nipmp 



(36) 
(37) 



where nip - 4-7TppSp/3 and equation dZTl i have been 
used. A large erosion rate is obtained when /j^ax is 
given by the planetesimal layer (requiring hp <«; hf) and 
when qp » 1. Under these conditions fragmentation by 
fragment-planetesimal collisions outcompetes fragmentation 
by planetesimal-planetesimal collisions. 

Adding the contribution from planetesimal-planetesimal 
collisions and fragment-planetesimal collisions, equations 
( |32] | and ( |36] |. gives a combined fragmentation rate of: 



-PF 



I + qp Tcoi.pp 



■qp 



hp_ 

T^col.PP \heff 



(38) 



2.5.4. Change in the characteristic size of the planetesimals and 
fragments. 

The planetesimal and fragment population are assigned a 
characteristic size {sf and sp). In the toy model it is possi- 
ble to follow the coagulation of planetesimals, increasing sp 
with time. If this feature is implemented (PG=1), planetesi- 
mals grow on a collision timescale (eq. Il26l ): 



dsp 
dt 



3T 



(PG=1), 



(39) 



col.PP 



Thus, for planetesimal collisions we can have the somewhat 
counterintuitive situation that they grow (increasing sp) while 
at the same time lose surface density due to fragmentation 

(eq. my 

Since bodies in the gravity regime become weaker when 
they are smaller, fragmentation of planetesimals triggers a 
collisional cascade. In our toy model the fragment mass (sf) 
represents the lower range of the cascade. As the initial value 



for i/r we choose the point where the strength and gravity con- 
tributions to Q*^ are equal (see Fig.|5]l. Thus, mass liberated 
at s = Sp flows instantaneously towards s = Sf. Typical ini- 
tial sizes for Sf are of the order of 100 m (for basalt or icy 
planetesimals). We only decrease Sf when fragment-fragment 
collisions become disruptive. The fragment size is then given 
by the implicit equation qpisf) = v^f/^Q*j{sf) - 0.5, where 
Uf is the random motion among fragments and the 0.5 v alue 
follows from comparison against lKobavashi et al.l (120101) . 

2.6. Treatment of solid's radial motions 

In protoplanetary disks, solids are expected to move radi- 
ally. Embryos may migrate due to type-I migration, planetes- 
imals can diffuse or scatter due to gravitational encounters, 
and radial orbital decays due to drag forces. In our context, 
radial drift of small particles (fragments) is especially impor- 
tant since it operates on short timescales. To capture this ef- 
fects, we will assign a characteristic drift timescale for the 
planetesimal and a fragment components, T^yk- In addition, 
the toy model can account for influx (accretion) from external 
regions (Mgxt + 0). Nevertheless, since our model is local, 
it only crudely captures the nature of these radial motions; a 
multidimensional extension (see below) is needed for a self- 
consistent treatment. 

Let us denote by Wsim the radial width of the local space 
around aq that is modeled by the toy model. The mass in 
solids, planetesimals {k - P) or fragments {k - F), is denoted 
Mw./t and changes according to 



; — = 2;raZ,, — ; — 

dt dt 



+ Mext.A- + Mint - 



Tdi-k 



(40) 



where the terms on the RHS of equation (l40b denote, respec- 
tively: (i) the increase in mass due to the growth of the feeding 
zone (the surface density outside Wsim may be different from 
E and is denoted I,'); (ii) the increase due to accretion from 
external regions; (iii) the change in mass within Wsim due to 
transfer of mass to a component other than k, e.g., fragmen- 
tation (eq. 1381 ) and embryo growth; and (iv) the mass loss 
from Wsim due to radial drift. Since the surface density is 
'^k = Miot,k/27TaWsim equation ( l40b translates into 



dLk_ 
dt 

with 



1 dW,„ 

Wsi^ dt 



-iK-^k) + 



ext,k 



InaW.i, 



r,k 



Vdi.k 



(RD=1) 



(41) 



(42) 



When radial drift is not implemented (RD=0), one can for- 
mally write Tdi.k = c»- Equation (l4Tl i very generally describes 
the evolution of 2/.. But the length scale Wsim yet needs to be 
specified. Here, we consider two choices for Wsim reflecting 
the Umits of a local or a global application of equation ( 1411 1. 

2.6.1. Local interpretation 

In a local model Wsim flo although Wsim should be cho- 
sen larger than the feeding/stirring zone of the embryo, i.e., 
Wsim > bRh. For such small Wsim, the drift timescale Tdr,k 
will generally be short (especially for fragments). Bound- 
ary effects, that is the mass flow Mext.i- coming into Wsim, are 
important and must be quantified. Ideally, Mext is obtained 
from the neighboring annulus as part of a multi-zone setup, 
Mext,*: - 2to'S'l'^^, where primes denote quantities from the 



Understanding how planets become massive 



11 



annulus exterior to Wsim- Although a multi-zone extension 
allows for a self-consistent treatment, it also considerably in- 
creases the overhead of the calculations. We will not imple- 
ment it here, but may choose do so in an upcoming work. 

Generally, therefore, we must specify the accretion rate, 
Me^t,k{t), as function of time. The picture here is of planet- 
formation going on in zone Wsim, with solids exterior to Wsim 
steadily drifting in. In § 14. 21 we consider such a scenario as a 
simple application of the toy model in its local setting. 

2.6.2. Global interpretation 

In a global setting, we assume that the physical condi- 
tions are similar for scales on the order of the disk radius, 
Wsim ~ flo- We no longer bother about mass in- or outflow due 
to growth of protoplanets or accretion from external regions 
(these are boundary effects). The first two terms of equation 
(HTl i are put 0. Specifically, we write 



dt 



TdY.k 
Qrift 



(43) 
(44) 



where Cddft ~ 1 is a fudge factor. We found that Cdrift - 
0.5 gives the best correspondence with multi-zone simulations 
models. 

The global extrapolation is tailored towards the oligarchic 
growth phase, since we can assume that neighboring oligarchs 
are of similar mass (although this assumption will break down 
on scales ~ao)- Thus, when the switch OL is set, equations 
( |43] | and ( l44l i apply. On the other hand, the local scenario (eqs. 
BTI and ll42l ) applies when we consider a single protoplanet 
(OL=0). 

2.6.3. The surface density evolution timescale 

Equation (l38T l is one of the terms contributing to the inter- 
nal mass flow lint in equations (HTt and ( |43] |: it re-arranges 
the mass among the three components, but conserves the total 
surface density in Wsim- Embryo growth also transfers mass; 
^kE - -C'^cc^E/Tgr,k- Note the reduction by C^^^ as, when 
oligarchy is considered, a third of the embryo growth comes 
from embryo-embryo collisions (see § 12.2. lb . 

Together, these two terms (i.e., eq. Il38l and E^e) give the 
total internal mass flow Eint,*. Let us for simplicity denote the 
first two terms of equation flTt by Sext.i, which represent the 
mass flow from exterior regions or from the expansion of the 
feeding zone. This term is therefore zero when 0L=1. Then, 
we can combine equations (HTt and (|43] | into 



■-ext,/; 



Tdr,k 



(45) 



where Td, ^ is given either by equation ( l42l i orl44l The surface 
density then evolves on a timescale of 



AnLk 



+ 1, 



2i 



1 

Tdi.k 



(46) 



Very often, when Tdr,k is short (i.e., for small fragments), a 
semi-steady state is reached where Ej; ^ (Eint,^ + Eext,*)7"dr,/t- 
In that case, \Tz,k\ » T^dr,* and it is computationally advanta- 
geous to evolve the evolution on timescales of Tz,k instead of 
the much shorter Td^k- 



2.7. Summary of the algorithm 

In essence, the algorithm consist of a loop in which the 
key quantities (the embryo mass Me, surface densities E^, and 
characteristic sizes st) are advanced over a timestep of arbi- 
trary length, until a final time (Tend) is reached. Within each 
iteration / we calculate: (i) the eccentricity of the fragment 
and size distribution at f = f,; (ii) the appropriate timestep Af,-; 
and (iii) the updated values for Me, "^k, and at r,+i - f,- + Af,-: 

1. Given an embryo mass Me, a particle size si^ and sur- 
face density Ej;, we calculate the growth mode - e qui- 
libri um o r non-equilibrium growth - as outlined in § 12.31 
and 12.41 This provides the embryo growth time Tgrk, 
the (Hill) eccentricity of the planetesimals and frag- 
ments, ehji, and the friction times Tfy^k- Using the cal- 
culated eccentricities, we obtain the fragmentation rates 
as outlined in § l2.5l and the timescales T^^k on which the 
fragment and planetesimal populations evolve, equa- 
tion ( |46] |. 

2. The next step concerns the choice of the timestep Af,. 
It is computationally efficient to choose a Af, as large 
as possible. Yet, we must resolve the system on the 
shortest of the timescales characterizing the evolution 
of the system: 

Af/ = ^ min {rgr,p; T^^y, T^y, T^y, T'coI.pp} , (47) 

where ^ is a control parameter (we choose ^ = 0.1 
throughout). Note that Tdr,F does not enter equation 
( |47] |. as it typically is very short. Instead, we resolve 
the evolution of the surface density. The planetesimal- 
planetesimal collision timescaleTcoi.pp is included to 
resolve planetesimal growth (§ 12. 5.4b . 

3. After the timestep is assigned, the quantities are ad- 
vanced and f,+i = f, + Af,. The embryo mass is increased 
by an amount AMe = MEAti/Tgr where Tgr is the total 
growth timescale due to accretion of planetesimals and 
fragments, T^/ = T^^'p + T^g/p- The planetesimal and 
fragment sizes ar e chan ged according to the prescrip- 
tion outlined in § 12.5.41 Finally, we update the surface 
densities of embryos, planetesimals, and fragments, us- 
ing Af,- and T^.i. 



3. COMPARISON AGAINST PREVIOUS STUDIES 

We test the toy model against previous literature studies. 
Throughout this and following section, we use the values 
listed in the second column of Table |2] for the parameters of 
the toy mode l, unless specified otherw ise. In particularly, 
we adopt the iBenz & Asphaud ( 119991) material parameters 
for ices (the Qos, Qog, bs and bg parameters). Furthermore, 
we adopt the minimum mass surface nebula (MMSN) model 
(Weidenschilling 1977b: Hayashi et al. 1985) as our bench- 
mark for the surface density and temperature structure of the 
disk: 



E,,(fl) =4.0 X 10^ g cm-2 

2^ini(fl)=/g"d' 2^9(0); 

W=125k(— ) (-) ; 



1.5 



(48) 
(49) 

(50) 
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Table 3 

Switch and parameter values 



Section Figure 



Switcli value 



Parameters deviating from standard value or range 











AE 


ER 


FR 


ND 


OL 


PG 


RD TD 


TS 


3.1 


Figure 


6 
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X 


X 
















3.2 


Figure 
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X 







1 


1 





3.3 


Figure 


8 
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1 


1 







1 


1 





4.1 


Figure 




I 





1 


1 







1 


1 X 


X 


4.1 


Figure 
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X 




1 


1 X 
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Figure 
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1 
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1 
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Figure 
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1 


X 








1 X 






ao ■ 

Q'ss 



= 7 g cm - 

3.2 AU 

3-35 AU; 
= 10"* -0.1 



^ini = 1-10 km; ip ii 



1-100 km 



gos = gog = 
«hw = 2 - 54 m s"' 
sp = 10-' - 10'' cm; 



S/. = 0; Mi„i = 10"^ Me; Afext > 



Note. — Overview of the simulations in §[3]and|4]listing the behavior of the switches: indicates the switch is turned off; 1 that it is switched on; 
and 'x' that its influence is examined. See Table[T]for the abbreviation and description of the switches and Table|2]for the default parameter values. 
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Figure 6. (left) Sample runs at 5 AU, testing the influence of atmospheres and fragmentation on the growth of the embryo. The equilibrium solution for the 
dynamical state of the planetesimal population has been assumed and radial drift is switched off. The solid black line conesponds to the default run without 
fragmentation and atmosphere enhancement. The dashed curve corresponds to a run including planetesimal fragmentation but without accounting for erosive 
collisions (fragments colliding with planetesimals). The dotted curve includes the erosive term, (right) The surface densities in planetesimals, fragments, and 
embryos for the fragmentation run without erosion (black curves) and including erosion (gray curves). 



where /gd is the gas-solid mass ratio, assumed to equal 57 
beyond the snowline (where ices contribute to the surface 
density), /mmsn an enhancement factor for disks more mas- 
sive than the MMSN, and the stellar luminosity. Note that 
equations (I48]l-(l50ll are merely a disguise for the local nature 
of the model, i.e., the results depend on the local value of E, Q, 
Tg, only. From the MMSN model we calculate the isothermal 

sound speed - y/ksTg/iu, the disk scaleheight, Hg = Cj/Q, 

and the gas density at the midplane, pg - "Lg/ y/lKHg. 

The values of the switches and parameters deviating form 
their defaults are listed in Table |3] for each of the model 
runs. Except indicated otherwise, we adopt an embryo mass 
of Mini.E = 10"*' Me, assuming that by this stage the 2- 
component assumption characterizing oligarchy has materi- 
alized (see § [Tland im i. Although in reality the start of the 
oligarchi c phase depends on disk radius and initia l planetesi- 
mal size (llda & Makind[T99llOrmel et al.ll20T0ah . we do not 
expect that changes in Mini ^ will matter greatly for the out- 
come of the results. 

3.1. Comparison aeainst \ChamberA i200^} 

iChambersI (I2006L 12008) presents a model to follow the 
growth of embryos. In Chambers (2006) the basic model is in- 
troduced, whose elements are similar to ours: a 3-component 
model of planetesimals, fragments, and embryos. Slowly , 
the complexity of the model is increased: IChamberi (l2006h 
already includes the full calculation of the dynamical state, 
atmosphere enhancement of the capture radius, and radial 



drift. IChambersI (120081) includes a more detailed fragmenta- 
tion model (with multiple size bins) and adds type-1 migra- 
tion, thereby including a spatial dimension (the disk radius) 
to the model. 

We test our toy model against the 2006 (local) ver- 
sion of Chambers' model. This implies that we do not 
treat vertical diffusion of fragments, turbulent stirring, neb- 
ular drag effects, planetesimal growth, or erosive collisions: 
TD=TS=ND=ER=PG=RD=0 (see Table[3j. The model is con- 
ducted for a disk radius of 5 AU at a surface density 
of 7 gem-'. Furthermore, we only consider the equilib- 
rium solution (the viscous stirring timescale balances friction 
timescale) to obtain the dynamical state of planetesimals. Re- 
sults are presented in Fig.|6]and should be compared to Figs. 
6 and 8 of|Chambers (2006). 

With these assumptions our setups are nearly identical and 
we should expect to obtain very similar results. The dashed 
curves labeled '(default)' in Fig. [6^ shows the embryo mass as 
function of time for the case without planetesimal fragmenta- 
tion and without atmosphere enhancement. The dotted curve 
shows the results when atmospheres are included (AE=1). The 
run without atmospheres reaches the isolation mass Mjso only 
after f = 10^ yr, whereas the run that includes the increase 
of the collisional accretion rates due to atmosphe res growth 
accelerates after Me ~ O.IM^, in agreement with (IChambersI 
120061) . The isolation mass is the total solid mass of an embryo 
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Time [yr] Time [yr] 

Figure 7. (left) Mass of protoplanetary embiyos vs. time for a ran neglecting fragmentation (black curve ) and a run including fragmentation (gray curve). The 
dashed pa rts of th e curves indicate that planetesimal growth follows the non-equilibrium solution (§ l2.4.2t : whereas the solid parts correspond to the equilibrium 
solution (§ 12. 4. It . (right) The surface density in planetesimals (Y.p, thick solid curves), fragments (Sf , thick dashed curve), and the characteristic sizes (sp and 
if ; thin solid curves) for the run without fragmentation (black curves) and including fragmentation (gray curves). 



within an annulus of bRi,: 



3/2 



V3M; 



(51) 



which equals Miso = 11.8 Me for the parameters in Fig.|6] 
When the embryo reaches the isolation mass, its surface den- 
sity dominates, » S/r -H Y.p, rendering the oligarchic ap- 
proximations (i.e., ce - 0) inapplicable (see § 12. 11 1. But this 
only concerns the final doubling time in Me- 

In IChamber"i (l2006h fragments are assumed to lie in a dy- 
namically cold {ep — Q) thin layer. To mimic this effect, we 
fix the initial fragment size at 1 cm. Gas drag will then en- 
sure that the eccentricities are indeed negligible. The a ccre- 
tion efficiency of these particles is very large (see § 12.2. 3b . and 
because fragments are not assumed to drift (RD=0), fragmen- 
tation accelerates the growth significantly. The solid black 
curve in Fig.|6h shows this behavior: after t - 10^ yr (roughly 
the collision timescale of the planetesimals) growth rates in- 
crease dramatically and the isolation mass is reached very 
quickly. Figure |6j3, where the surface densities of the compo- 
nents are plotted further illustrates this point. Although frag- 
mentation is important, Y.p never becomes very large; frag- 
ments are very quickly accreted by the embryos. 

The gray solid curve in Fig.|6t and the gray curves in 
Fig.|6j3 correspond to the outcome of a run that includes the 
erosion term (i.e., the second term on the RHS of eq. Il38l ). 
This will be used in all of the following runs. As can be seen 
from Fig. the surface density in planetesimals now declines 
precipitously after t - 10^ yr. A positive feedback effect 
emerges: the higher the fragments' E/r, the higher cfLp/dt. In 
this laminar case (ass = 0), the sandblasting of planetesimals 
occurs very fast, because the cm-size fragments are confined 
to a thin layer at the disk's midplane. 



3.2. Comparison asainst lkobavashi et al.\ ( 120701) 

iKobavashi et al.l (120101) presents a, multi-zone, multi-bin, 
model to follow the core formation. They start from a 
monodisperse population of planetesimals and then compute 
the runaway growth and oligarchic growth stages. Their 
model is global: the total surface density of particles changes 
by radial drift motions from the outer disk to the inner 
disks. Planetesimal fragmentation accelerates this inward 
flow of solids. Indeed, iKobavashi et al.. (.2010.) found that 



fragmentation was generally not conducive to growth, a 
somewhat diff'erent con clusion from Chambers (2006, 2008). 
iKobayashi et al.l (120101) only considers laminar disks, and 
also ignore nebular drag effects. Thus, the ID , TS and ND 
switches are all turned off". The PG switch is turned on, 
whereas the FR switch can be turned on or off, depending on 
the model test. IKobayashi et al.l (1201 ll) include the effect of 
atmospheres. 

We focus on the standard model of lKobayashi et al.l (120101) . 
where the MMSN is assumed. Figure shows the evolu- 
tion of the embryo mass at a dis k radius of 3.2 AU. Thi s plot 
should be compared to Fig. 7 of lKobavashi et al] ( 120101) . The 
initial, dashed part of the curves indicate that planetesimals 
follow the non-equilibrium solution, whereas solid curves in- 
dicates the equilibrium solution (see § I2.4| |. The black curve 
is a run without fragmentation; the gray curve includes the 
effects of planetesimal fragmentation and erosion. 

The growth-only run (black curve) follows a similar pattern 
as the corresponding curve of lKobayashi et all (120101) . After 
3 X 10^ yr the growth mode switches from non-equilibrium 
growth to equilibrium growth. An increase in the growth rate 
can then be expected since for the equilibrium regime gravita- 
tional focusing factors no longer decrease. However, no steep 
increase is seen due to two opposing effects. First, planetes- 
imal growth increases the mass of the planetesimals, which 
renders gas friction less effective and results in a smaller fo- 
cusing factor (or larger e;,). Second, towards the end of the 
simulation, the planetesimal mass reservoir becomes empty. 
Figure]?}) illustrates these effects: sp (thin solid, black line) 
increases and "Lp (thick solid, black line) decreases towards 
the end of the simulation. 

A somewhat more complex behavior can be seen when 
fragmentation is switched on (FR=1; gray curves in Fig.|7^,b). 
On timescales of ~Tco\ pp planetesimals collide and fragment. 
The size of the f ragments, determined from the criterion out- 
lined in § 12.5.41 is initially ~100 m and further decreases as 
fragments collide among themselves. When their surface den- 
sity is large, accretion of fragments can become very effective 
since their eccentricities are low. The surface density in frag- 
ments peaks towards f = 10^ yr at a; 10% of the initial sur- 
face density in fragments. At this point the accretion rate is 
substantial. However, the fragmentation and decreasing sp 
cause a stronger orbital decay as Tf,- = Orf, decreases towards 
unity (see eq. ETI ). As a result, I,p quickly decreases. In 
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Figure 8. The embryos' masses after lO' yr as function of disk radius for models that include fragmentation. Open black circles coirespond to models including 
atmosphere enhancement of the collisional cross section (AE=1), while open black squares correspond to runs without atmospheres. Note that each of these 
points corresponds to a separate run of our local model. The soli d gray symbols give the corresponding results of the multi-zone simulations offKobayashi et J] 
1201 II) . The solid line indicates the isolation mass M^so (eq- 1511 ). The nine panels each vary according to the disk mass with respect to the MMSN (eq. 1481 ) 
and the initial size of the planetesimals .sp jni . 



the end, then, the combined effect of fragmentation and radial 
drift stalls the growth as it depletes the solid mass reservoir. 

3.3. Including Atmospheres (AE-1) 

We next run the toy model at several disk radii, pretending, 
to mimic the setting of a global model. We will test the effects 
of the initial size of planetesimals (ip.ini), the inclusion of at- 
mospheres (the AE switch), and the disk mass (which trans- 
lates in a larger gas density pg and initial solid surface density 
Sini). In total 144 models are run. In Fig. [8] we plot the final 
mass (after 10^ yr) of the runs. Black open squares corre- 
sponds to runs without atmosphere enhancement and black 
open circles to runs with AE=1. These simulations should be 
compared to Figs. 3-5 of Kobayashi et al. (201 1), which data 
are shown in Fig.[8]by solid gray symbols. 

Although the general correspondence is good, occasion- 
ally the_predictionsof_fli^ are seen to deviate from 
the iKobavashi et al.l (120101 1201 Ih simulations. The embryo 
masses for the 3xMMSN, 1 km runs at 10- 20 AU in the toy 
model, for example, are lower than those of iKobayashi et al] 
( |201 lii) . but the discrepancy especially becomes severe after 
AE is switched on. This is understandable, since, as we saw in 
§ 13.11 growth can accelerate dramatically when AE becomes 
important. The large discrepancy therefore merely reflects the 
sensitivity to the onset of this effect - rather than indicating a 
fundamental flaw in the model. 



A more systematic difference is that for the toy model em- 
bryo masses show a stronger dependence on the disk mas (en- 
hancement over MMSN). Whereas the final embryo masses 
of the toy model mostly lag those of Kobayashi et al. (1201 lb 
at IxM MSN, the toy model has overtaken the lKobavashi et alj 
( 1201 ll) simulations for the lOxMMSN runs, especially when 
AE is switched on. This discrepancy may partly reflect 
our simplistic treatment of the atmosphere models (see 
Appendix|B]), where we assume that the opacity stay constant 
(and low) - assumptions that may break down for massive 
atmospheres. More fundamentally, strong drift motions of 
solids causes the breakdown of the oligarchy approximation: 
i.e., that most of the mass is in planetesimals and fragments, 
rather than embryos. For example, in the lOxMMSN, 1 km run 
planetesimal fragmentation is very efficient and the embryos 
wil l dominate the surface d ensity sooner rather than later. In 
the lKobavashi et all ( 1201 lb simulations the embryos then ex- 
cite themselves to large eccentricities^ which frustrates their 
growth. In the toy model, however, eg = is enforced and 
embryos always merge while keeping a separation of bR/,- 

Altogether, the match of the toy model to th e more so- 
phisticated simulations by IKobavashi et all ( 1201 lb is very sat- 

The simulations conducted by IKobavashi et alj 1201 11) did not in- 
clude eccentric ity d amping due to gravitational torques with the gas disk 
iTanaka et al. 2003) . 
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Figure 9. Sensitivity of embryogrowth to the gas and material properties and ND,TS,TD switches. In each panel the thick gray line con'esponds to the 5AU, 
3xMMSN, 1 km run from Fig.[8]with ER=FR=OL=PG=RD= 1 and AE=ND=TS=TD=8. In (a) we switch on turbulence by activating the TS and TD switches, and 
examine the sensitivity of the results to the turbulence-ffss parameter, whose standard value is ffss = IQ-'*. In (foj we put 2oj = Qog = 0, indicative of a rabble-pile 
nature of the planetesimals without any internal strength. Collisions produce a copious amount of mm-size fragments, whose accretion behavior is very sensitive 
to the AE, ND, and TD switches. In (c) we reduce the disk headwind o^w from its standard value of 54 m s-'. 



isfactory. The general trends are similar, the final masses 
agree within a factor of 3 for most runs, and large discrepan- 
cies are mainly caused by the boost embryos experience due 
to crossing a certain mass threshold for the atmosphere ef- 
fects kick in. Our resul ts reconfirm the main conclusion from 
iKobavashi et alj (1201 lb : that the planetesimal fragmentation- 
radial drift tandem makes it hard to grow large cores, even 
when atmosphere treatment is accounted for. A great advan- 
tag e of the toy model is th at it is very fast: whereas it takes 
the IKobavashi et all (1201 ll) simulations 5-10 CPU hours on a 
modern PC to complete a single disk model, the eight corre- 
sponding runs of the toy model are finished in »10 seconds. 

4. APPLICATIONS 

Due to its flexibility, there are countless imaginable applica- 
tions that can be envisioned with the toy model. A full inves- 
tigation is beyond the scope of this paper. Here, we illustrate 
the toy model by two examples: (i) a modest sensitivity study 
to investigate in what way the results of the previous section 
depend on adopted gas and material properties; and (ii) the 
calculation of the accretion efficiency of fragments, drifting 
past a single embryo. In both examples, we focus primarily 
on the behavior of small particles, as determined by the AE, 
ND, and TD switches. 

4.1. Sensitivity to gas and material properties 

Here we take a particular run of Fig. [8]- 5AU, 3xMMSN, 10 
km, no atmosphere enhancement - and explore how its results 
change if we change the adopted gas and material properties: 
the level of turbulence in the disk (Fig.|9^); the material prop- 
erties (Fig. |9j') and the disk headwind parameter ^hw (Fig. |9]:). 
In each panel the reference model of Fig.[8]is indicated by the 
thick gray line. 

In the toy model, turbulence is controlled by the dimension- 
less ass parameter and is activated by setting one of the TS or 
TD switches. Turbulent stirring (TS) due to density inhomo- 
geneities in the gas may outweigh over viscous stirring when 
ffss is large, or when embryos are small (see AppendixICTl. 
Turbulent diffusion (TD) causes the particle layer to puff up, 
thereby decreasing the local densities of particles near the 
midplane. Figure|9] shows that for a turbulence parameter 
of QTss - iO '^ (solid curve) growth is initially very slow. 
The slow growth is due to turbulent stirring, which excites 
the planetesimals to large, but constant eccentricities. Over 



time, the Hill eccentricity e/, decreases as the embryo grows; 
growth speeds up after a slow start and viscous stirring takes 
over eventually. This 'delayed' onset of (runaway) growth has 
been seen before ( lOrmel et al.l20105lGressel et alj20Tll) and 
depends on the severity of the stirring, here parametrized by 

In addition, large ass will cause planetesimals to fragment 
'prematurely', i.e., when the embryo is still small. Fragments 
drift away before the embryo has an opportunity to accrete 
them, resulting in much smaller embryos. This effect is re- 
flected in the ass = 10"^^ (dashed curve) and the ass = 10"^ 
(where embryo growth does not take off) runs. There is a cru- 
cial difference between fragmentation triggered by turbulent 
stirring and viscous stirring. When viscous stirring is impor- 
tant, the embryo is by definition massive and it will always 
sweep up some fraction of the solids. But turbulent stirring 
is insensitive to the embryo mass, which renders it potentially 
much more harmful to embryo growth. 

The size of the fragments and their production rate are 
largely determined by the material properties. In Fig.|9l5 we 
consider the case that the planetesimals are rubble piles by 
putting Qos - Q()a = 0, leaving the gravitational binding en- 
ergy as the sole source of strength (eq. 1291 ; Fig.|5]l. Further- 
more, the fragment size is fixed at 1 mm at all times. This 
setup is very conducive to growth, because the radial drift for 
mm-size particles is not so effective as it is for m-size parti- 
cles. As a result, embryos in the standard model (solid line) 
accrete these particles before they drift away; and embryos 
grow large on a very short timescale (~10 yr!). Nebular 
drag effects (ND=1) cause the fragments to be dragged along 
with the gas, rendering it somewhat more difficult to accrete 
them, especially when the protoplanet is small (dashed curve). 
However, when the gas is only slightly turbulent (ass = 10""*; 
dotted curve) accretion virtually shuts down since the small 
particles are now distributed over a large scaleheight: embryo 
growth stalls at masses <0.1 M®. When small particles con- 
tribute significantly to the mass in solids, core formation will 
be very sensitive to the state of the turbulence. 

Finally, we vary the value of the disk headwind Lih„. This 
parameter primarily determines the drift timescale, and hence 
the residence time of particles near the embryo. Figure|9]: 
shows that a reduction by a factor of 5 (from the default value 
of Dhw = 54 m s"' to 10 m s"') results in embryos that are a 
factor of three more massive. Reducing the headwind param- 
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eter even further reinforces these findings. The 100 m-size 
'fragments' stay longer in the accretion zone, which enhances 
the likelihood that the embryo can capture them. A larger 
capture radius due to atmospheres (AE=1) exacerbates these 
effects (dotted curve in Fig.|9};). Although in Fig.|9}; frag- 
ments are typically ~10 m in size, the reduction of Hhw favors 
embryo growth for all fragment sizes. 

The value of the disk headwind Uhw, or, equivalently, the 
pressure gradient (eq. ifTTI ). can be affected in several ways: 

1. Variations of the global pressure profile, i.e., of the sur- 
face density and temperature index. From equation (fTTb 
one can show that Uhw depends on the exponents of the 
temperature and surface density of the gas. Thus, devi- 
ations from the MMSN value p = -3/2 (as in X oc a'') 
will be directly reflected in The same holds for a 
different temperature structure, e.g., in optically thick 
disks. 

2. Local changes in dP/da. More potent changes arise lo- 
cally, if we account for variations in t he accretion rat e 
of the gas, triggered by dead zones (lGammielll996h . 
The interfaces between MRl dead and active zones in 
particular are promising candidates in this respect. At 
these interfaces, the kinematic viscosity vj- of the tur- 
bulence changes abruptly. To maintain a constant gas 
accretion rate [Mg cc vj^g) a decrease with disk ra- 
dius in vt is accompanied by an increase in I,g and a 
corresponding increase in the local pressure gradient 
(Kretke & Lin 2007). Possibly, this even reverses the 
sign of the pressure gradient, such that drifting parti- 
cles will accumulate in high pressure maxima. 

3. Collective effects. When dust particles are sufficiently 
concentrated they will act collectively and drag the gas 
with them in the direction of Keplerian velocities, in ef- 
fect reducing the headwind a dust particle feels. Here, 
we have not accounted for collective effect in our ex- 
pressions for the radial drift (eqs. ||2T| and ll22l ). b ut this 
is straightforward to i mplement ([Ta naka et al.j l2005t 
lEstrada & Cuzzil 120081: iBai & Stonell2010a.) . Particle- 
hydrodynamic si mulations often displ ay strong collec- 
tive effects (e.g.. lJohansen et all2009t) . 

4.2. Efficiency of accreting inward-drifting particles by a 
single embryo 

State of the art particle-hydrodyna mical simulation s show 
that planetesima ls may form big (iJohansen et al.l l2009t 
ICuzzi et aI1l2010l) . perhaps with radii up to 10^ km, out of a 
'sea' of small part icles - an idea that also finds support in the 
asteroidal record (Bottke et alj 120051: iMorbidelU etall 120091: 
but see Weidenschilling 2011] for an alternative interpreta- 
tion). In this context, it is very relevant to s tudy the interac- 
tion between embryos and particles directly. IXie et al.) (1201 Oh 
studied the sweeping-up of small particles by km-size plan- 
etesimals, and showed that the planetesimals acquired mass 
pretty fast, even though gravitational focusing was omitted. 
Here, we extend their study to the regime of embryos. 

We consider a scenario where particles from the outer disk 
drift past a single protoplanet situated at 5 AU. For simplic- 
ity, we do not treat a planetesimal component (Zp - 0); 
the planetesimal-related switches (FR, ER, PG and TS) are 
therefore irrelevant. The solids that drift inwards are all small 
particles of fixed size Sf. As we consider a single protoplanet 
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Figure 10. Accretion efficiency of inward-drifting particles by a single 
embryo. Each symbol corresponds to a run where a total mass of 100 in 
fragments of size Sf (x-axis) drifts past an embryo over a timescale of ~10* 
yr. The fragment size does not evolve. Curves coiTespond to choices of the 
AE, ND, and TD switches: AE=ND=TD=9 (black solid line), AE=ND=1 (dashed 
solid line), AE=ND=TD=1 (gray lines with ffss = 10"^ and a^s = 10"'). 

the OL switch is off and the surface density evolution is de- 
scribed by equation ( |4T] ). The accretion rate Mext is given by 

Mex,(f) = — exp[-f/f,], (52) 

tc 

where Mtot is the total, time-integrated, amount of solids 
that drifts by and tc the characteristic timescale. Here, we 
choose Mtot - 100 and t^ = 10^ yr. It must be said 
that these parameters and equation ( [52] i are completely ad- 
ho c and purely chosen for i llustrative convenience. We refer 
to lYoudin & Chiand (120041) and Birnstiel et al. (subm.) for 
physically-motivated prescriptions. The initial mass of the 
embryo is fixed at Mi^i E = 0.01 M^. 

The embryo mass after f » f^h provides us with the ac- 
cretion efficiency of the embryo. This quantity is plotted in 
Fig.[TO]as function of fragment size. Each symbol in this plot 
denotes a different run with the symbols connected by lines 
to guide the eye. The default suite of models (diamonds) rep- 
resents the case without accounting for drag effects or tur- 
bulent diffusion (AE=ND=TD=0). One ontices that big boul- 
ders of 10-100 m in size have axl% (time-averaged) proba- 
bility of being accreted. Somewhat smaller, m-size particles 
fare much worse, though. They are optimally coupled to the 
gas (Tfr = 1) and speed by at the maximum drift velocity of 
fdr ~ fhw = 54 m s"' (see eq. 1211 ). Smaller particles, how- 
ever, reside longer in the feeding zone of the planet and have 
a much higher probability of getting accreted - for mm-size 
particles it becomes almost 100%. 

When we allow for drag-effects (i.e., an atmosphere- 
enhanced capture radius, and nebular drag; dashed curve in 
Fig.fTOl). the accretion rate increases across the spectrum of 
sizes. Heavier fragments now reach accretion probabilities of 
several tens of percents, sufficient to form large cores. Accre- 
tion efficiencies also remain large for the smallest sizes, due 
to their long residence times in combination with their low 
scaleheights. Turbulence, which will stir up these small par- 
ticles, will therefore tend to negate these effects, see the gray 
curves in Fig.[TOl For a (modest) a^s value of 10""* it already 
shows much lower accretion efficiency for mm and cm-size 
particles; for a^s = 10"^ the accretion efficiency is even lower 
These finding reflects those of Fig.|9}5. 

In conclusion. Fig. [TOl shows that the accretion efficiency 
of fragments can become sufficiently large to build sizable 
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cores {xlQ M^), although for smaller particles this depends 
heavily on the turbulence levels. MRJ dead zones may there- 
fore favor core growth, although recent s tudies indicate that 
particles are nonetheless very diffusive (iTurner et alJ 120 lOt 
lOkuzumi & Hiros^l201 lb . In addition, the results will depend 
on the headwind velocity, i.e., the (local) pressure gradient 
of the gas. FigurefTOl corresponds to Vhvi - 54 ms"'; but, 
as we saw from Fi g. [9}:, a lower headwi nd will enhance the 
accretion rate. The lJohansen et all (1201 ll) simulations, for ex- 
amples, show that m-size boulders concentrate heavily. These 
'particle bands' are dense enough to drag the gas along, such 
that the headwind vanishes. Collective effect will mitigate low 
efficiencies seen in Fig.lTOlfor Tfr = 1 particles. 

5. SUMMARY AND OUTLOOK 

We have presented a new model to follow the protoplanet 
growth stage in planet formation. The model is a toy model: 
it assumes that the mass distribution can be approximated by 
three components (embryos, planetesimals, and fragments) 
and does not involve a radial dimension. To obtain the equa- 
tions that govern the evolution, we have employed simple 
physically-based recipes. Yet the toy model is intended to 
be general: a wide array of physical relevant processes can be 
included at the discretion of the user. In particular, we have 
focused on a more accurate treatment of fragments, for which 
coupling to gas is important. Due to our modular setup (the 
switches) the toy model is intended to expose the key mecha- 
nisms that drive protoplanet growth. 

One should be aware of the limitations, as we have out- 
lined in § 12.11 The (at most) three component approximations 
entails that the toy model cannot model the runaway growth 
phase (the phase preceding oligarchy) and does a bad job in 
modeling the fragmentation cascade. The oligarchic assump- 
tions will also break down when the mass in embryos (Zg) 
starts to exceed that of planetesimals and fragments - a sit- 
uation that may be reached sooner rather than later in cases 
where drift is important. Furthermore, as of now, type I mi- 
gration, or planetesima l driven migration jKirsh et al.ll2009t : 
ICapobianco et al.l201 IL Ormel & Ida in prep) are not included 
either. Our model is in it strictest sense only locally valid. 

We have shown in § [3] that our r e sults are in good agree- 
ment with those of lKobavashi et al.l (1201011201 11) . These cal- 
culations indicate that planetesimal fragmentation (a natu- 
ral outcome of oligarchic growth) is quite harmful to proto- 
planet growth, due to their strong orbital decay. This con- 
clusion contrasts with many earlier works, which suggest 
that small planetesimals, or fragments, are very conducive 
for growth (Rafikov 2004; Goldreich et al. 200 4 iChambersI 
I200d iKenvon & Bromlev..2009; .Guilera et al.u201 Ih . How- 
ever, Chambersl(l2008h alreadv warned that fragmentation acts 
like a double-edged sword - increasing collision rates but at 
the same time depleting the disk of so lids due to orbital de- 
cay - whereas iKobavashi et al.l (1201 1[) finds that fragmenta- 
tion is detrimental to embryo growth. For this reason, they 
prefer bigger planetesimals since these are stronger and less 
collisional - both effects reduce the fragmentation rate. But 
the disadvantage of big planetesimals is that they are damped 
poorly and growth is relatively slow. 

But it may still be premature to provide a final answer 
to the question whether small planetesimals, or small frag- 
ments, promote or hinder (fast) growth. Very small parti- 
cles have the additional handicap that turbulence stirs them 
up as we saw in Fig.|9j5 and Fig.[TO] Few works concern- 
ing protoplanet growth have accounted for turbulent stirring 



dBromlev & Kenvonll201 ll are a recent exception), which is 
somewhat pecuhar since qtss is a well known - and criti- 
cal - param eter for the pre-planetesim al growth state (e.g., 
iBirnstiel et a l. 2010; Zsom et al.ir201 lb . Indeed, fr om an ob- 
servation perspective mm/cm particles dominate (the mass) 
of protoplanetary disks and are the likely precursors of the 
first generation of big plan e tesimals e.g. prot oplanet seeds 
dJohansen et al.1 12007l l2009t ICuzzi et al.l 120101) . Therefore, 
we better investigate the accretion behavior of smal l particles 
(IJohansen & Lacerdal l20Tot lOrmel & Klahd l20To[ Fig.\M- 
Here we reemphasize that our treatment of atmosphere en- 
hancements (AE) and nebular drag effects (ND) is somewhat 
'ad-hoc' (see our comments towards the end of § 12.2.41 ). 

In a follow-up study we intend to return to the question re- 
garding the efficiency of accretion of (big) planetesimal vs. 
(small) fragment. The result of this work offers some tenta- 
tive clues, though. First, the efficiency scales with the mag- 
nitude of the headwind ^hw - the deviation of the orbital ve- 
locity of the gas from Keplerian. It seems very hard to grow 
big protoplanets (~10 M®) out of small particles when Hhw 
is not reduced or with some form of replenishment of solids. 
Interestingly, a similar conclusion may hold for planetesimal 
formation itself (Cuzzi et al. 2010; Bai & Stone 2010b). 

A second important point is that the accretion efficiency 
is strongly size dependent (Fig.fTOb. (Sub)km-size planetes- 
imals accrete well, but are prone to collisional fragmentation. 
The resulting m-size fragments fare much worse though in 
terms of accretion efficiency, despite the fact that they have a 
large specific accretion rate. Strong radial drift for these parti- 
cles is the limiting factor. The accretion efficiency of smaller 
particles strongly depends on the turbulent state of the disk. 
Overall, these sensitivities indicate that accretion is very much 
dependent on the condition regarding the gaseous state (pres- 
sure gradient, turbulence) and on the material strength of the 
particles, which to a large extent determines the typical frag- 
ment size. 

In the end, any growth scenario must be tested against ob- 
servati onal constraints, e.g., t h e size distribution of t he as- 
teroid dMorbidelh et all 120091: IWeidenschili in3l2011b dust 
production rates in deb r is dis ks ( Shann on & Wul 1201 It 
[Kenvon & Bromlev" '2008', 20 10|), o r the exoplanet census 
dAlibertet al. 2011; Ida & Linii2010l) . For detailed compari- 
son, however, our toy model with its discrete, three compo- 
nent approximation is a rather crude tool. But it is a great 
tool to conduct a preliminary exploration of the vast parame- 
ter space. This provides the user with a first, crude compar- 
ison with the observational data and a means to identify the 
relevant physical processes that are at work, thereby guiding 
further sophisticated modeling efforts. We believe that such 
an interplay between a toy- and sophisticated model is a very 
powerful asset. 
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APPENDIX 

A. PcoL IN THE DRAG-MODIFIED REGIME (ND=1) 

For completeness, we provide a summary of the algorithm that lOrmel & KlahU (120101) (OKIO) provide to derive Pcoi when the 
equations of motion of the (small) particles are significantly affected by gas drag. For this work we have taken the opportunity to 
streamline the fitting expressions for Pcoh slightly improving the correspondence to the numerical integrationsQ 

In the drag-modified regime, Pcoi is a function of three (dimensionless) parameters: the dimensionless friction time rfr - Tf^Q.; 
the 'headwind parameter' (u: = fhw/^h (which is via i>/, a function of planet mass); and the planet size oe = Re/Ri,- The 
OKIO expressions cover three regimes, where the gas drag-gravity interaction is qualitatively different. From small to large 
these are: (i) the full settling regime; (ii) the hyperbolic regime; (iii) and the three body regime. The regime boundaries are 
Tfr = T* - min(12/f^„ 2) and Tfr = fu,. The hyperbolic regime disappears for large protoplanets. 

For each of these regimes OKIO gives expressions for the impact radius b^oi and the approach velocity Va- For the settling 

' The modifications with respect to OKIO entail the definition of the ex- 
pression T* and the extension of an exponentially-declining tail to bj^. 
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Table 4 

Expressions for Pcol when ND=1. 



Regime 


Condition 


Impact radius" focol 


Impact velocity «„ 


Smoothed impact radius, 


Setthng (set) 


Tf, < T* 


Equation lAH 


3h^J2 + 4, 


fo,e,exp[-(rfr/r')''«] 


Hyperbohc (hyp) 


r* < Tfr < 






^hyp 


Three body (3b) 


Tf,. > max[rfr, f„] 


Uoe + 1.0/Tf,. 


3.2 


fejb exp[-(0.7f,/Tfr)5] 



Note. — Summary of expressions to derive Pcol = ^fafe'^i when nebular drag influences the particle trajectories around 
the protoplanet (ND=1). Expressions and parameters are normalized with Hill units: Tfr = r^n, i^^ = «hw/wA, and qe = 
Re/R/i- The critical friction times is defined by r* = min(12/f[„, 2). 
" The minimum impact radius is always the geometrical radius, bcai = 0!e- 



regime, valid for Tfr < t* one obtains bco\ by solving a cubic equation 

+ ^b'- - 8Tfr = 0, (Al) 

whose real solution is denoted feset- Similarly, we find the impact radii for the hyperbolic and three-body regime (^hyp and /73b, 
respectively). Each of these regimes, furthermore, has a different expression for the approach velocity Vg. These expressions are 
summarized in Table |4] 

The expressions for ZjcoI and Va have been derived by OKIO using a simple physical model. These expressions are not continuous 
across boundaries. In particular, this pertains to the behavior of fehyp which is often much lower than bsei and /73b near the 
boundaries. To provide for a continuous transition one can smooth the impact radii fecoi by appending ^set and /73b with an 
exponentially-decaying tail, see the last column of TableU] This causes the regimes to overlap, e.g., the validity of the full 
settling regime is extended beyond > tJ^.. The smoothing exponentials are empirically determined by fitting the resulting Pcoi 
to the Pcoi obtained from numerical integration. 

The 2D accretion rate is then given by Pcoi = 2^*01"" where b*^^ is the empirical (smoothed) variant of the analytically-derived 
bco\ (see TablelU. Thus, one obtains Pset - '^Ket^^'et and, similarly, f hyp and P3b. By virtue of the exponential continuation, 
multiple solutions now exist for a given rfr and for the final Pcoi we take the maximum: 

Peol = max (P,,e„ Phyp, Pih) ■ (A2) 
B. CALCULATION OF THE RADIUS ENHANCEMENT BY ATMOSPHERES, ASSUMING CONSTANT OPACITY 

We will assume that the putative atmosphere of a growing protoplanet is radiatively supported. The equations for stellar 
structure then read: 

P=-^pr (Bl) 
dP GM 

— --— P (B2) 



dr r 
dT 3kLc p 



dr ^AncTsB r^T^ 



(B3) 



where P is the pressure, kg is Boltzmann's constant, //mn the molecular weight, p the (gas) density, T the temperature, r the 
height of the atmosphere as measured from the center of the protoplanet, k the opacity, Lc the luminosity of the protoplanet, ctsb 
Stefan-Boltzmann's constant and G Newton's gravitational constant. We normalize the radial coordinate r to the Bondi radius, 
X - r/Ri,, where Rb = GMIyc], - Pglpg the isothermal sound speed of the nebula, and j - 1.4. Equa tions (IBll i-( IB3l l are 
further normalized by the nebula quantities, i.e., p = P/Pg, cr - pipg, and 9 - T/Tg. Equations (IBlb - (IB3t , in nondimensional 
form, then read 



p = cr9 (B4) 
dp cr 
dx 



where Wneb is the only remaining parameter 3 



g 

^ The definition of W^ej, dift'ers by a factor of 4 from the Wo defined by 
llnaba & Ik oma lOO^,. 
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Figure 11. Solutions for the density cr = p/pg us function of dimensionless radius x in case of constant opacity throughout the atmosphere. At large x the density 
follows the isothermal curve. H owever, after a critical radius xi (indicated by the cross) the envelope becomes pressure-supported and the density follows a x^^ 
law. Dashed-gray curves are the llnaba & Ikomj 120031) solution for the density structure. 

Dividing equation ( IB6l l by equation iB5i we can readily solve for the temperature as function of the pressure 

04 = 1 + 4W„ebO^ - 1), (B8) 

where the nebular boundary condition 9(p = 1) = 1 has been applied. We now approximate equation dBSI l assuming is small 
(0 <K 1; nearly isothermal) or large (0 » 1; pressure dominated). The boundary between these regimes is set at a radius xi and a 
corresponding density cri. By construction, both solutions join at this point. 

B.l. Limit AW nehP 1 (nearly isothermal) 
When 4WnebP ^ 1 equation ( |B8t can be approximated as 6* a: 1 + W^ehip ~ !)■ We then obtain 



P^CtQ^ (r{\ + WnebO^ - 1)) * H^nebtr' + ^^(1 - W„eb). 

With this equation and the chain rule, we write equation ( IB5l l as 



Integration gives 



dcr cr 

ax x^ 



2W„,b(r + (1 - Wneb) \og(r^^ + C, 



Here the integration constant Cg, obtained from cr(x = 1) = 1, evaluates to Cg = 2Wneb - 7- We therefore find 



1 ^ ^ ^ 2Wneb(tr- l) + logCr 
X J 



(B9) 

(BIO) 
(BID 

(B12) 



(where we assumed Wneb ^ !)■ 



B.2. Limit 4WnehP ^ 1 (pressure dominated) 

In this case we approximate equation (IB8b as ff* x AWmhP- Then, equation (IB6b reduces to dO/dx = -yWnebp/0^x^ - -y/4-x^ 
and 6 - ylAx + C\, where C\ is another integration constant. Using equation (IB4b we have; 

^ 1 / 



P 

(T - — 



4W„eb 4W„eb 

The integration constant Ci may be found from the condition o-(xi) = o"!, i.e., Ci + y/A-xi = (4WnebO'i)'^"' and 

cr 



^'neb \ Ax Axxj 



4W„eb 



(B14) 



Thus, for X «; min(l,xi) the density scales as the cube of inverse radius, in agreement with previous studies (IStevensonll 19821 : 
ITiiaba & Ikoma.2003.) . 

B.3. Results 

We yet need to specify the transition between the nearly-isothermal and the pressure-domi nate d regimes, that is a-\ (and 
corresponding p\). If we put the transition at 4Wneb/5i = 1, cri can be obtained from equation (IB9b : cr\ ^ /?i/(l + WnehP\) ~ 
We therefore define 



5W„eb 



(B15) 
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Figure 12. Solutions for the capture radius normalized to the core radius, as function of particle size for several protoplanet masses. The protoplanet is 
placed at a disk radius of 5 AU and accretes particles at a rate of 1 Me Myr"' . The eccentricity of the particles is fixed at e/, = 4, independent of their size. 



(from which x\ follow s from eg. IIB12|). Th e density structure calculated by our method is in excellent agreement with the full 
analytic solution from (Inab_a_&lkoma 2003). In Fig.[TT]we show several examples. 

Our approximation provides a 1-1 relation between the radius and the density, cr - cr(x). Ilnaba & Ikomal (120031) calculate that 
a particle must experience a peak density of 

(6 + el)si,ps 

Pa = ' (B16) 

in order to lose a sufficient amount of its (3-body) energy to become captured by the protoplanet. Inverting the expressions 
for (t(x) then gives us a direct solution for the capture radius (xa) in closed form. That is, using equation ( IB12b for densities 
o-a = Pa/Pg < o"i and 

l.l + l(4W„eb)'^^K/3-^n (B17) 

for densities cr^, larger than (j\ . 

In Fig.[T2]we plot the radius enhancement factor of the embryo (Rg/Rc) as function of particle size protoplanets of mass 0.1, 
1.0, and 10 M®. For simplicity, we have fixed the Hill eccentricity at 4e/,, although in reality this will be a function of particle 
radius too. Likewise, the accretion rate is fixed at M = 1 yr ' . With these parameters the Wneb values are 2.8 x 10"'*, 1 .3 x lO"'*, 
and 5.9 x 10"^, respectively. The enhancement or Ra is largest for smaller particles as they are most affected by the drag. The 
increase in Ra continues until the 'boundary' of the atmosphere is hit, which is gi ven by the Bondi radius. The radius increase is 
also a steep function of protoplanet mass. Figure[T2]can be compared to Fig. 2 of lKobavashi et al.l ( 1201 lb . 

C. STIRRING BY TURBULENT DENSITY FLUCTUATIONS 

Planets (and planetesimals) are scattered by gas (over)densities induced by turbulence. The statistical distribution of these 
torque fluctuations is determined by two key quantities: the amplitude of the rms-fluctuations, crp and the correlation time Tc- 
These allow us to define a diffusion coefficient, Dj - ct^Tc such that the induced eccentricity change becomes 

Aj jD]i 

— ^ (CI) 
J J 

where j = a^Q is the specific angular momentum. Squaring and differentiating with respect to time gives the stirring rate 

The magnitude of the fluctuations, crp, is redefined in terms of a nondimensionless parameter y,: 

CI.y,a'n^ 

(Tr = — , (C3) 

dBamteau & Linll2010l) with C = 2.4 x 10~. If we assume that = O ' we find 



de^ ICy,a^I., 



dt \ Mi, 



Q.. (C4) 
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iBaruteau & LinI ( 120101) relate y, to the diffusion parameter ass, 7i - 8.5 x lO^-alpHg/a. Inserting this expression in equation 
iC4i . these numerical constants we obtain a stirring rate of 

-4xl0-ass -7:r^ " (C5) 



dt \ 

and a turbulent stirring timescale of: 



2^2 25xio-3///,fli,r' , 



It is instructive to compare the turbulent and viscous stirring timescales (eq. 



5 X 10-^ e'-R„ I Hga,i:A -^ ^ 10^^ Me i Hgapllg ^^^^ 



/ ass r' Pvs / Mg/M. y//g/fl\-V WM. '^-^ 

"^•°lTo^j ToI^o^Utt] 1^0^) ' ^^^^ 

where we used that e/e;, = Rh/ao - (M^/SM*)^. Thus, turbulent stirring may over dominate viscous stirring for small protoplan- 
ets (small Me), massive disks (large I.g) at large disk radii (where the flaring is stronger), and large ass- 



